GPU accelerated investigation of a dual-frequency driven nonlinear oscillator
Ferenc Hegedűs^{1}, Werner Lauterborn^{2}, Ulrich Parlitz^{3} and Robert Mettin^{2}
^{1}Department of Hydrodynamic Systems, Budapest University of Technology and Economics, Budapest, Hungary
^{2}Drittes Physikalisches Institut, Georg-August-Universität Göttingen, Göttingen, Germany
^{3}Research Group Biomedical Physics, Max Planck Institute for Dynamics and Self-Organization, Göttingen, Germany and Institute for Nonlinear Dynamics, Georg-August-Universität Göttingen, Göttingen, Germany
Summary. The bifurcation structure of a dual-frequency driven, second order nonlinear oscillator (Keller–Miksis equation) is investigated by exploiting the high computational resources of professional GPUs. The numerical scheme of the applied initial value problem solver was the explicit, adaptive Runge–Kutta–Cash–Karp method with embedded error estimation using solutions of order 4 and 5. The four dimensional parameter space (amplitudes and frequencies of the driving) is explored by means of several high resolution bi-parametric plots with the amplitudes as control parameters at fixed frequencies. The resolution of the control parameter plane is 500 × 500 with 10 initial conditions at each parameter pair (altogether 2.5 million initial value problems in each bi-parametric plot). The program code for one fine parameter scan runs approximately 50 times faster on a Tesla K20 GPU (Kepler architecture) than on an Intel i7-4790 4 core CPU even applying double precision floating point operations.
High resolution bifurcation structure
The bifurcation structure of nonlinear harmonic oscillators has been intensively studied during the last decades. As a function of a single control parameter, the dynamics of such systems is relatively well understood and the accumulated knowledge is summarized in many textbooks [1]. The real challenge today is to extend the numerical investigation to multi-dimensional parameter space. This requires orders of magnitude larger computational resources to determine the fine bifurcation substructure of nonlinear oscillators confidently. Therefore, many numerical studies on high resolution multi-dimensional parameter scans investigate iterated maps [2] or continuous systems where e.g. a simple 4th order Runge–Kutta scheme with fixed time step is sufficient [3]. In our case, however, an adaptive initial value problem solver is mandatory due to the possible very different time scales of the solutions of our model. This makes the computations more resource intensive. In order to obtain results within reasonable time, the very high floating point processing power of professional GPUs were exploited.
The employed nonlinear oscillator was the Keller–Miksis equation
which is a second order ordinary differential equation describing the radial oscillation of a single spherical bubble, for details see [4]. Here 𝑅(𝑡) is the time dependent bubble radius. The parameter values during the computations (gas bubble in water) were as follows: liquid density 𝜌_{𝐿} = 997.1 kg/m^{3}, sound speed 𝑐_{𝐿} = 1497.3 m/s and viscosity 𝜇_{𝐿} = 8.902^{−4} Pa s; vapour pressure 𝑝_{𝑉} = 3166.8 Pa; surface tension 𝜎 = 0.072 N/m; ambient (static) pressure 𝑃_{∞} = 1 bar; bubble size 𝑅_{𝐸} = 10 μm; polytropic exponent 𝑛 = 1.4 (adiabatic behaviour). The dual-frequency driving of the system
𝑝_{∞}(𝑡)=𝑝_{𝐴1}sin(𝜔_{1}𝑡)+𝑝_{𝐴2}sin (𝜔_{2}𝑡)
is a time varying pressure field, where 𝑝_{𝐴1} and 𝑝_{𝐴2} are the pressure amplitudes (also control parameters); 𝜔_{1} and 𝜔_{2} are the corresponding angular frequencies. During the simulations, the frequencies were normalized by the linear undamped resonance frequency 𝜔_{0} = 340 kHz of the system [4] (𝜔_{𝑟𝑒𝑙} = 𝜔/𝜔_{0}). The model was solved by an adaptive and explicit Runge–Kutta–Cash–Karp initial value problem solver with error estimation of orders 4 and 5. To obtain reliable data for a thorough topological analysis of the bifurcation structure of the system, several high resolution bi-parametric plots were simulated in the 𝑝_{𝐴1}−𝑝_{𝐴2} bi-parametric plane with 10 randomly chosen initial conditions at each parameter pair to reveal co-existing attractors. The resolution of the pressure amplitudes was Δ𝑝 = 0.01 bar in both directions (500 × 500 grid points). After the initial transient, the properties of a found attractor were recorded (Poincaré points, Lyapunov exponent, winding number etc.). Due to the explicit numerical scheme (only function evaluations are needed), the parameter studies of our problem (solution of 2.5 million independent initial value problems for each bi-parametric plot) are well suited for GPUs. The applied relative frequencies were selected from 𝜔𝑟𝑒𝑙 = 1/10, 1/5, 1/3, 1/2, 1, 2, 3, 5, 10. Simulations were performed at every possible combinations of the frequency values, which means altogether 36 high resolution bi-parametric plots. Results at some pairs of relative frequencies are shown in Fig. 1, where the Lyapunov exponent is plotted as a function of the control parameters 𝑝_{𝐴1} and 𝑝_{𝐴2}. The periodic (negative) and chaotic (positive Lyapunov exponent) solutions are indicated by the greyscale and coloured domains, respectively. Although the presented diagrams are valuable for further investigation, it must be emphasized that the present study focuses only on performing the numerical simulations rather that the topological description of the plotted bifurcation structure itself. The available GPUs were a GeForce Titan Black (Kepler architecture), 2 Tesla K20 (Kepler architecture) and 5 Tesla M2050 (Fermi architecture). The simulation times were approximately 50 times faster on a Kepler architecture and 25 times faster on a Fermi architecture than on an Intel i7-4790 4 core CPU, using double precision floating point operations.
Figure 1: High resolution Lyapunov exponent diagrams in the 𝑝𝐴1−𝑝𝐴2 parameter plane describing the topology of simple periodic (greyscale) and chaotic (colour) domains at different relative frequency combinations.
Acknowledgement
This research was supported by the Deutsche Forschungsgemeinschaft (DFG) grant no. ME 1645/7-1, and by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences.
References
[1] Strogatz, S. H. (2014) Nonlinear Dynamics and Chaos with Applications to Physics, Biology, Chemistry, and Engineering. Westview Press, Boulder, Colorado.
[2] de Souza, S. L. T. and Lima, A. A. and Caldas, I. L. and Medrano-T, R. O. and Guimaraes-Filho, Z. O. (2012) Self-similarities of periodic structures for a discrete model of a two-gene system. Phys. Lett. A 376:1290-1294.
[3] Gallas, J. A. C. (2015) Periodic oscillations of the forced Brusselator. Mod. Phys. Lett. B 29:1530018.
[4] Lauterborn, W. and Kurz, T. (2010) Physics of bubble oscillations. Rep. Prog. Phys. 73:106501.