1 /*
2 * chaotic_system.cpp
3 *
4 * This example demonstrates how one can use odeint to determine the Lyapunov
5 * exponents of a chaotic system namely the well known Lorenz system. Furthermore,
6 * it shows how odeint interacts with boost.range.
7 *
8 * Copyright 2011-2012 Karsten Ahnert
9 * Copyright 2011-2013 Mario Mulansky
10 *
11 * Distributed under the Boost Software License, Version 1.0.
12 * (See accompanying file LICENSE_1_0.txt or
13 * copy at http://www.boost.org/LICENSE_1_0.txt)
14 */
15
16
17 #include <iostream>
18 #include <boost/array.hpp>
19
20 #include <boost/numeric/odeint.hpp>
21
22 #include "gram_schmidt.hpp"
23
24 using namespace std;
25 using namespace boost::numeric::odeint;
26
27
28 const double sigma = 10.0;
29 const double R = 28.0;
30 const double b = 8.0 / 3.0;
31
32 //[ system_function_without_perturbations
33 struct lorenz
34 {
35 template< class State , class Deriv >
operator ()lorenz36 void operator()( const State &x_ , Deriv &dxdt_ , double t ) const
37 {
38 typename boost::range_iterator< const State >::type x = boost::begin( x_ );
39 typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
40
41 dxdt[0] = sigma * ( x[1] - x[0] );
42 dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
43 dxdt[2] = -b * x[2] + x[0] * x[1];
44 }
45 };
46 //]
47
48
49
50 //[ system_function_with_perturbations
51 const size_t n = 3;
52 const size_t num_of_lyap = 3;
53 const size_t N = n + n*num_of_lyap;
54
55 typedef boost::array< double , N > state_type;
56 typedef boost::array< double , num_of_lyap > lyap_type;
57
lorenz_with_lyap(const state_type & x,state_type & dxdt,double t)58 void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
59 {
60 lorenz()( x , dxdt , t );
61
62 for( size_t l=0 ; l<num_of_lyap ; ++l )
63 {
64 const double *pert = x.begin() + 3 + l * 3;
65 double *dpert = dxdt.begin() + 3 + l * 3;
66 dpert[0] = - sigma * pert[0] + 10.0 * pert[1];
67 dpert[1] = ( R - x[2] ) * pert[0] - pert[1] - x[0] * pert[2];
68 dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2];
69 }
70 }
71 //]
72
73
74
75
76
main(int argc,char ** argv)77 int main( int argc , char **argv )
78 {
79 state_type x;
80 lyap_type lyap;
81
82 fill( x.begin() , x.end() , 0.0 );
83 x[0] = 10.0 ; x[1] = 10.0 ; x[2] = 5.0;
84
85 const double dt = 0.01;
86
87 //[ integrate_transients_with_range
88 // explicitly choose range_algebra to override default choice of array_algebra
89 runge_kutta4< state_type , double , state_type , double , range_algebra > rk4;
90
91 // perform 10000 transient steps
92 integrate_n_steps( rk4 , lorenz() , std::make_pair( x.begin() , x.begin() + n ) , 0.0 , dt , 10000 );
93 //]
94
95 //[ lyapunov_full_code
96 fill( x.begin()+n , x.end() , 0.0 );
97 for( size_t i=0 ; i<num_of_lyap ; ++i ) x[n+n*i+i] = 1.0;
98 fill( lyap.begin() , lyap.end() , 0.0 );
99
100 double t = 0.0;
101 size_t count = 0;
102 while( true )
103 {
104
105 t = integrate_n_steps( rk4 , lorenz_with_lyap , x , t , dt , 100 );
106 gram_schmidt< num_of_lyap >( x , lyap , n );
107 ++count;
108
109 if( !(count % 100000) )
110 {
111 cout << t;
112 for( size_t i=0 ; i<num_of_lyap ; ++i ) cout << "\t" << lyap[i] / t ;
113 cout << endl;
114 }
115 }
116 //]
117
118 return 0;
119 }
120