1 /*
2 * elliptic_functions.cpp
3 *
4 * Copyright 2011-2013 Mario Mulansky
5 * Copyright 2011-2012 Karsten Ahnert
6 *
7 * Distributed under the Boost Software License, Version 1.0.
8 * (See accompanying file LICENSE_1_0.txt or
9 * copy at http://www.boost.org/LICENSE_1_0.txt)
10 */
11
12
13
14
15 #include <iostream>
16 #include <fstream>
17 #include <cmath>
18
19 #include <boost/array.hpp>
20
21 #include <boost/numeric/odeint/config.hpp>
22
23 #include <boost/numeric/odeint.hpp>
24 #include <boost/numeric/odeint/stepper/bulirsch_stoer.hpp>
25 #include <boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp>
26
27 using namespace std;
28 using namespace boost::numeric::odeint;
29
30 typedef boost::array< double , 3 > state_type;
31
32 /*
33 * x1' = x2*x3
34 * x2' = -x1*x3
35 * x3' = -m*x1*x2
36 */
37
rhs(const state_type & x,state_type & dxdt,const double t)38 void rhs( const state_type &x , state_type &dxdt , const double t )
39 {
40 static const double m = 0.51;
41
42 dxdt[0] = x[1]*x[2];
43 dxdt[1] = -x[0]*x[2];
44 dxdt[2] = -m*x[0]*x[1];
45 }
46
47 ofstream out;
48
write_out(const state_type & x,const double t)49 void write_out( const state_type &x , const double t )
50 {
51 out << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] << endl;
52 }
53
main()54 int main()
55 {
56 bulirsch_stoer_dense_out< state_type > stepper( 1E-9 , 1E-9 , 1.0 , 0.0 );
57
58 state_type x1 = {{ 0.0 , 1.0 , 1.0 }};
59
60 double t = 0.0;
61 double dt = 0.01;
62
63 out.open( "elliptic1.dat" );
64 out.precision(16);
65 integrate_const( stepper , rhs , x1 , t , 100.0 , dt , write_out );
66 out.close();
67
68 state_type x2 = {{ 0.0 , 1.0 , 1.0 }};
69
70 out.open( "elliptic2.dat" );
71 out.precision(16);
72 integrate_adaptive( stepper , rhs , x2 , t , 100.0 , dt , write_out );
73 out.close();
74
75 typedef runge_kutta_dopri5< state_type > dopri5_type;
76 typedef controlled_runge_kutta< dopri5_type > controlled_dopri5_type;
77 typedef dense_output_runge_kutta< controlled_dopri5_type > dense_output_dopri5_type;
78
79
80 dense_output_dopri5_type dopri5 = make_dense_output( 1E-9 , 1E-9 , dopri5_type() );
81 //dense_output_dopri5_type dopri5( controlled_dopri5_type( default_error_checker< double >( 1E-9 , 0.0 , 0.0 , 0.0 ) ) );
82
83 state_type x3 = {{ 0.0 , 1.0 , 1.0 }};
84
85 out.open( "elliptic3.dat" );
86 out.precision(16);
87 integrate_adaptive( dopri5 , rhs , x3 , t , 100.0 , dt , write_out );
88 out.close();
89 }
90