1 /*
2 * Copyright 2011-2013 Mario Mulansky
3 * Copyright 2012 Karsten Ahnert
4 *
5 * Distributed under the Boost Software License, Version 1.0.
6 * (See accompanying file LICENSE_1_0.txt or
7 * copy at http://www.boost.org/LICENSE_1_0.txt)
8 *
9 * Example for the lorenz system with a 3D point type
10 */
11
12 #include <iostream>
13 #include <cmath>
14
15 #include <boost/operators.hpp>
16
17 #include <boost/numeric/odeint.hpp>
18
19
20 //[point3D
21 class point3D :
22 boost::additive1< point3D ,
23 boost::additive2< point3D , double ,
24 boost::multiplicative2< point3D , double > > >
25 {
26 public:
27
28 double x , y , z;
29
point3D()30 point3D()
31 : x( 0.0 ) , y( 0.0 ) , z( 0.0 )
32 { }
33
point3D(const double val)34 point3D( const double val )
35 : x( val ) , y( val ) , z( val )
36 { }
37
point3D(const double _x,const double _y,const double _z)38 point3D( const double _x , const double _y , const double _z )
39 : x( _x ) , y( _y ) , z( _z )
40 { }
41
operator +=(const point3D & p)42 point3D& operator+=( const point3D &p )
43 {
44 x += p.x; y += p.y; z += p.z;
45 return *this;
46 }
47
operator *=(const double a)48 point3D& operator*=( const double a )
49 {
50 x *= a; y *= a; z *= a;
51 return *this;
52 }
53
54 };
55 //]
56
57 //[point3D_abs_div
58 // only required for steppers with error control
operator /(const point3D & p1,const point3D & p2)59 point3D operator/( const point3D &p1 , const point3D &p2 )
60 {
61 return point3D( p1.x/p2.x , p1.y/p2.y , p1.z/p2.z );
62 }
63
abs(const point3D & p)64 point3D abs( const point3D &p )
65 {
66 return point3D( std::abs(p.x) , std::abs(p.y) , std::abs(p.z) );
67 }
68 //]
69
70 //[point3D_norm
71 // also only for steppers with error control
72 namespace boost { namespace numeric { namespace odeint {
73 template<>
74 struct vector_space_norm_inf< point3D >
75 {
76 typedef double result_type;
operator ()boost::numeric::odeint::vector_space_norm_inf77 double operator()( const point3D &p ) const
78 {
79 using std::max;
80 using std::abs;
81 return max( max( abs( p.x ) , abs( p.y ) ) , abs( p.z ) );
82 }
83 };
84 } } }
85 //]
86
operator <<(std::ostream & out,const point3D & p)87 std::ostream& operator<<( std::ostream &out , const point3D &p )
88 {
89 out << p.x << " " << p.y << " " << p.z;
90 return out;
91 }
92
93 //[point3D_main
94 const double sigma = 10.0;
95 const double R = 28.0;
96 const double b = 8.0 / 3.0;
97
lorenz(const point3D & x,point3D & dxdt,const double t)98 void lorenz( const point3D &x , point3D &dxdt , const double t )
99 {
100 dxdt.x = sigma * ( x.y - x.x );
101 dxdt.y = R * x.x - x.y - x.x * x.z;
102 dxdt.z = -b * x.z + x.x * x.y;
103 }
104
105 using namespace boost::numeric::odeint;
106
main()107 int main()
108 {
109
110 point3D x( 10.0 , 5.0 , 5.0 );
111 // point type defines it's own operators -> use vector_space_algebra !
112 typedef runge_kutta_dopri5< point3D , double , point3D ,
113 double , vector_space_algebra > stepper;
114 int steps = integrate_adaptive( make_controlled<stepper>( 1E-10 , 1E-10 ) , lorenz , x ,
115 0.0 , 10.0 , 0.1 );
116 std::cout << x << std::endl;
117 std::cout << "steps: " << steps << std::endl;
118 }
119 //]
120