1[/============================================================================ 2 Boost.odeint 3 4 Copyright 2011-2012 Karsten Ahnert 5 Copyright 2011-2013 Mario Mulansky 6 Copyright 2012 Sylwester Arabas 7 8 Use, modification and distribution is subject to the Boost Software License, 9 Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at 10 http://www.boost.org/LICENSE_1_0.txt) 11=============================================================================/] 12 13 14[section Chaotic systems and Lyapunov exponents] 15 16[import ../examples/chaotic_system.cpp] 17 18In this example we present application of odeint to investigation of the properties of chaotic 19deterministic systems. In mathematical terms chaotic refers to an exponential 20growth of perturbations ['__delta x]. In order to observe this exponential growth one usually solves the equations for the tangential dynamics which is again an ordinary differential equation. These equations are linear but time dependent and can be obtained via 21 22['d __delta x / dt = J(x) __delta x] 23 24where ['J] is the Jacobian of the system under consideration. ['__delta x] can 25also be interpreted as a perturbation of the original system. In principle 26['n] of these perturbations exist, they form a hypercube and evolve in the 27time. The Lyapunov exponents are then defined as logarithmic growth rates of 28the perturbations. If one Lyapunov exponent is larger then zero the nearby 29trajectories diverge exponentially hence they are chaotic. If the largest 30Lyapunov exponent is zero one is usually faced with periodic motion. In the 31case of a largest Lyapunov exponent smaller then zero convergence to a 32fixed point is expected. More information's about Lyapunov exponents and nonlinear 33dynamical systems can be found in many textbooks, see for example: E. Ott "Chaos is 34Dynamical Systems", Cambridge. 35 36To calculate the Lyapunov exponents numerically one usually solves the equations of motion for ['n] perturbations and orthonormalizes them every ['k] steps. The Lyapunov exponent is the average of the logarithm of the stretching factor of each perturbation. 37 38To demonstrate how one can use odeint to determine the Lyapunov exponents we choose the Lorenz system. It is one of the most studied dynamical systems in the nonlinear dynamics community. For the standard parameters it possesses a strange attractor with non-integer dimension. The Lyapunov exponents take values of approximately 0.9, 0 and -12. 39 40The implementation of the Lorenz system is 41 42`` 43const double sigma = 10.0; 44const double R = 28.0; 45const double b = 8.0 / 3.0; 46 47typedef boost::array< double , 3 > lorenz_state_type; 48 49void lorenz( const lorenz_state_type &x , lorenz_state_type &dxdt , double t ) 50{ 51 dxdt[0] = sigma * ( x[1] - x[0] ); 52 dxdt[1] = R * x[0] - x[1] - x[0] * x[2]; 53 dxdt[2] = -b * x[2] + x[0] * x[1]; 54} 55`` 56We need also to integrate the set of the perturbations. This is done in parallel to the original system, hence within one system function. Of course, we want to use the above definition of the Lorenz system, hence the definition of the system function including the Lorenz system itself and the perturbation could look like: 57 58`` 59const size_t n = 3; 60const size_t num_of_lyap = 3; 61const size_t N = n + n*num_of_lyap; 62 63typedef std::tr1::array< double , N > state_type; 64typedef std::tr1::array< double , num_of_lyap > lyap_type; 65 66void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t ) 67{ 68 lorenz( x , dxdt , t ); 69 70 for( size_t l=0 ; l<num_of_lyap ; ++l ) 71 { 72 const double *pert = x.begin() + 3 + l * 3; 73 double *dpert = dxdt.begin() + 3 + l * 3; 74 dpert[0] = - sigma * pert[0] + 10.0 * pert[1]; 75 dpert[1] = ( R - x[2] ) * pert[0] - pert[1] - x[0] * pert[2]; 76 dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2]; 77 } 78} 79`` 80 81The perturbations are stored linearly in the `state_type` behind the state of the Lorenz system. 82The problem of '''lorenz()''' and '''lorenz_with_lyap()''' having different state types may be solved putting the Lorenz system inside a functor with templatized arguments: 83 84`` 85struct lorenz 86{ 87 template< class StateIn , class StateOut , class Value > 88 void operator()( const StateIn &x , StateOut &dxdt , Value t ) 89 { 90 dxdt[0] = sigma * ( x[1] - x[0] ); 91 dxdt[1] = R * x[0] - x[1] - x[0] * x[2]; 92 dxdt[2] = -b * x[2] + x[0] * x[1]; 93 } 94}; 95 96void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t ) 97{ 98 lorenz()( x , dxdt , t ); 99 ... 100} 101 102`` 103This works fine and `lorenz_with_lyap` can be used for example via 104`` 105state_type x; 106// initialize x.. 107 108explicit_rk4< state_type > rk4; 109integrate_n_steps( rk4 , lorenz_with_lyap , x , 0.0 , 0.01 , 1000 ); 110`` 111This code snippet performs 1000 steps with constant step size 0.01. 112 113A real world use case for the calculation of the Lyapunov exponents of Lorenz system would always include some transient steps, just to ensure that the current state lies on the attractor, hence it would look like 114 115`` 116state_type x; 117// initialize x 118explicit_rk4< state_type > rk4; 119integrate_n_steps( rk4 , lorenz , x , 0.0 , 0.01 , 1000 ); 120`` 121The problem is now, that `x` is the full state containing also the 122perturbations and `integrate_n_steps` does not know that it should only use 3 123elements. In detail, odeint and its steppers determine the length of the 124system under consideration by determining the length of the state. In the 125classical solvers, e.g. from Numerical Recipes, the problem was solved by 126pointer to the state and an appropriate length, something similar to 127 128`` 129void lorenz( double* x , double *dxdt , double t, void* params ) 130{ 131 ... 132} 133 134int system_length = 3; 135rk4( x , system_length , t , dt , lorenz ); 136`` 137 138But odeint supports a similar and much more sophisticated concept: __boost_range. To make the steppers and the system ready to work with __boost_range the system has to be changed: 139 140[system_function_without_perturbations] 141 142This is in principle all. Now, we only have to call `integrate_n_steps` with a 143range including only the first 3 components of ['x]: 144 145[integrate_transients_with_range] 146 147[note Note that when using __boost_range, we have to explicitly configure the 148stepper to use the `range_algebra` as otherwise odeint would automatically 149chose the `array_algebra`, which is incompatible with the usage of __boost_range, because the original state_type is an `array`.] 150 151Having integrated a sufficient number of transients steps we are now able to calculate the Lyapunov exponents: 152 153# Initialize the perturbations. They are stored linearly behind the state of the Lorenz system. The perturbations are initialized such that [' p [subl ij] = __delta [subl ij]], where ['p [subl ij]] is the ['j]-component of the ['i].-th perturbation and ['__delta [subl ij]] is the Kronecker symbol. 154# Integrate 100 steps of the full system with perturbations 155# Orthonormalize the perturbation using Gram-Schmidt orthonormalization algorithm. 156# Repeat step 2 and 3. Every 10000 steps write the current Lyapunov exponent. 157 158[lyapunov_full_code] 159 160The full code can be found here: [github_link examples/chaotic_system.cpp chaotic_system.cpp] 161 162[endsect] 163