• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*
2  * stepper_details.cpp
3  *
4  * This example demonstrates some details about the steppers in odeint.
5  *
6  *
7  * Copyright 2011-2012 Karsten Ahnert
8  * Copyright 2012 Mario Mulansky
9  * Copyright 2013 Pascal Germroth
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 #include <iostream>
17 #include <boost/array.hpp>
18 #include <boost/bind.hpp>
19 #include <boost/numeric/odeint.hpp>
20 
21 using namespace std;
22 using namespace boost::numeric::odeint;
23 
24 const size_t N = 3;
25 
26 typedef boost::array< double , N > state_type;
27 
28 //[ system_function_structure
sys(const state_type &,state_type &,const double)29 void sys( const state_type & /*x*/ , state_type & /*dxdt*/ , const double /*t*/ )
30 {
31     // ...
32 }
33 //]
34 
sys1(const state_type &,state_type &,const double)35 void sys1( const state_type &/*x*/ , state_type &/*dxdt*/ , const double /*t*/ )
36 {
37 }
38 
sys2(const state_type &,state_type &,const double)39 void sys2( const state_type &/*x*/ , state_type &/*dxdt*/ , const double /*t*/ )
40 {
41 }
42 
43 
44 //[ symplectic_stepper_detail_system_function
45 typedef boost::array< double , 1 > vector_type;
46 
47 
48 struct harm_osc_f1
49 {
operator ()harm_osc_f150     void operator()( const vector_type &p , vector_type &dqdt )
51     {
52         dqdt[0] = p[0];
53     }
54 };
55 
56 struct harm_osc_f2
57 {
operator ()harm_osc_f258     void operator()( const vector_type &q , vector_type &dpdt )
59     {
60         dpdt[0] = -q[0];
61     }
62 };
63 //]
64 
65 //[ symplectic_stepper_detail_system_class
66 struct harm_osc
67 {
f1harm_osc68     void f1( const vector_type &p , vector_type &dqdt ) const
69     {
70         dqdt[0] = p[0];
71     }
72 
f2harm_osc73     void f2( const vector_type &q , vector_type &dpdt ) const
74     {
75         dpdt[0] = -q[0];
76     }
77 };
78 //]
79 
main(int argc,char ** argv)80 int main( int argc , char **argv )
81 {
82     using namespace std;
83 
84     // Explicit stepper example
85     {
86         double t( 0.0 ) , dt( 0.1 );
87         state_type in , out , dxdtin , inout;
88         //[ explicit_stepper_detail_example
89         runge_kutta4< state_type > rk;
90         rk.do_step( sys1 , inout , t , dt );               // In-place transformation of inout
91         rk.do_step( sys2 , inout , t , dt );               // call with different system: Ok
92         rk.do_step( sys1 , in , t , out , dt );            // Out-of-place transformation
93         rk.do_step( sys1 , inout , dxdtin , t , dt );      // In-place tranformation of inout
94         rk.do_step( sys1 , in , dxdtin , t , out , dt );   // Out-of-place transformation
95         //]
96     }
97 
98 
99 
100     // FSAL stepper example
101     {
102         double t( 0.0 ) , dt( 0.1 );
103         state_type in , in2 , in3 , out , dxdtin , dxdtout , inout , dxdtinout;
104         //[ fsal_stepper_detail_example
105         runge_kutta_dopri5< state_type > rk;
106         rk.do_step( sys1 , in , t , out , dt );
107         rk.do_step( sys2 , in , t , out , dt );         // DONT do this, sys1 is assumed
108 
109         rk.do_step( sys2 , in2 , t , out , dt );
110         rk.do_step( sys2 , in3 , t , out , dt );        // DONT do this, in2 is assumed
111 
112         rk.do_step( sys1 , inout , dxdtinout , t , dt );
113         rk.do_step( sys2 , inout , dxdtinout , t , dt );           // Ok, internal derivative is not used, dxdtinout is updated
114 
115         rk.do_step( sys1 , in , dxdtin , t , out , dxdtout , dt );
116         rk.do_step( sys2 , in , dxdtin , t , out , dxdtout , dt ); // Ok, internal derivative is not used
117         //]
118     }
119 
120 
121     // Symplectic harmonic oscillator example
122     {
123         double t( 0.0 ) , dt( 0.1 );
124         //[ symplectic_stepper_detail_example
125         pair< vector_type , vector_type > x;
126         x.first[0] = 1.0; x.second[0] = 0.0;
127         symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
128         rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , x , t , dt );
129         //]
130 
131         //[ symplectic_stepper_detail_system_class_example
132         harm_osc h;
133         rkn.do_step( make_pair( boost::bind( &harm_osc::f1 , h , _1 , _2 ) , boost::bind( &harm_osc::f2 , h , _1 , _2 ) ) ,
134                 x , t , dt );
135         //]
136     }
137 
138     // Simplified harmonic oscillator example
139     {
140         double t = 0.0, dt = 0.1;
141         //[ simplified_symplectic_stepper_example
142         pair< vector_type , vector_type > x;
143         x.first[0] = 1.0; x.second[0] = 0.0;
144         symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
145         rkn.do_step( harm_osc_f1() , x , t , dt );
146         //]
147 
148         vector_type q = {{ 1.0 }} , p = {{ 0.0 }};
149         //[ symplectic_stepper_detail_ref_usage
150         rkn.do_step( harm_osc_f1() , make_pair( boost::ref( q ) , boost::ref( p ) ) , t , dt );
151         rkn.do_step( harm_osc_f1() , q , p , t , dt );
152         rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , q , p , t , dt );
153         //]
154     }
155 
156     // adams_bashforth_moulton stepper example
157     {
158         double t = 0.0 , dt = 0.1;
159         state_type inout;
160         //[ multistep_detail_example
161         adams_bashforth_moulton< 5 , state_type > abm;
162         abm.initialize( sys , inout , t , dt );
163         abm.do_step( sys , inout , t , dt );
164         //]
165 
166         //[ multistep_detail_own_stepper_initialization
167         abm.initialize( runge_kutta_fehlberg78< state_type >() , sys , inout , t , dt );
168         //]
169     }
170 
171 
172 
173     // dense output stepper examples
174     {
175         double t = 0.0 , dt = 0.1;
176         state_type in;
177         //[ dense_output_detail_example
178         dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > > dense;
179         dense.initialize( in , t , dt );
180         pair< double , double > times = dense.do_step( sys );
181         (void)times;
182         //]
183 
184         state_type inout;
185         double t_start = 0.0 , t_end = 1.0;
186         //[ dense_output_detail_generation1
187         typedef boost::numeric::odeint::result_of::make_dense_output<
188             runge_kutta_dopri5< state_type > >::type dense_stepper_type;
189         dense_stepper_type dense2 = make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() );
190         (void)dense2;
191         //]
192 
193         //[ dense_output_detail_generation2
194         integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() ) , sys , inout , t_start , t_end , dt );
195         //]
196     }
197 
198 
199     return 0;
200 }
201