ceph/src/boost/libs/numeric/odeint/examples/lorenz_point.cpp

   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
  30     point3D()
  31         : x( 0.0 ) , y( 0.0 ) , z( 0.0 )
  32     { }
  33
  34     point3D( const double val )
  35         : x( val ) , y( val ) , z( val )
  36     { }
  37
  38     point3D( const double _x , const double _y , const double _z  )
  39         : x( _x ) , y( _y ) , z( _z )
  40     { }
  41
  42     point3D& operator+=( const point3D &p )
  43     {
  44         x += p.x; y += p.y; z += p.z;
  45         return *this;
  46     }
  47
  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
  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
  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;
  77     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
  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
  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
 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 //]