120 lines
2.7 KiB
C++
120 lines
2.7 KiB
C++
/*
|
|
* Copyright 2011-2013 Mario Mulansky
|
|
* Copyright 2012 Karsten Ahnert
|
|
*
|
|
* Distributed under the Boost Software License, Version 1.0.
|
|
* (See accompanying file LICENSE_1_0.txt or
|
|
* copy at http://www.boost.org/LICENSE_1_0.txt)
|
|
*
|
|
* Example for the lorenz system with a 3D point type
|
|
*/
|
|
|
|
#include <iostream>
|
|
#include <cmath>
|
|
|
|
#include <boost/operators.hpp>
|
|
|
|
#include <boost/numeric/odeint.hpp>
|
|
|
|
|
|
//[point3D
|
|
class point3D :
|
|
boost::additive1< point3D ,
|
|
boost::additive2< point3D , double ,
|
|
boost::multiplicative2< point3D , double > > >
|
|
{
|
|
public:
|
|
|
|
double x , y , z;
|
|
|
|
point3D()
|
|
: x( 0.0 ) , y( 0.0 ) , z( 0.0 )
|
|
{ }
|
|
|
|
point3D( const double val )
|
|
: x( val ) , y( val ) , z( val )
|
|
{ }
|
|
|
|
point3D( const double _x , const double _y , const double _z )
|
|
: x( _x ) , y( _y ) , z( _z )
|
|
{ }
|
|
|
|
point3D& operator+=( const point3D &p )
|
|
{
|
|
x += p.x; y += p.y; z += p.z;
|
|
return *this;
|
|
}
|
|
|
|
point3D& operator*=( const double a )
|
|
{
|
|
x *= a; y *= a; z *= a;
|
|
return *this;
|
|
}
|
|
|
|
};
|
|
//]
|
|
|
|
//[point3D_abs_div
|
|
// only required for steppers with error control
|
|
point3D operator/( const point3D &p1 , const point3D &p2 )
|
|
{
|
|
return point3D( p1.x/p2.x , p1.y/p2.y , p1.z/p2.z );
|
|
}
|
|
|
|
point3D abs( const point3D &p )
|
|
{
|
|
return point3D( std::abs(p.x) , std::abs(p.y) , std::abs(p.z) );
|
|
}
|
|
//]
|
|
|
|
//[point3D_norm
|
|
// also only for steppers with error control
|
|
namespace boost { namespace numeric { namespace odeint {
|
|
template<>
|
|
struct vector_space_norm_inf< point3D >
|
|
{
|
|
typedef double result_type;
|
|
double operator()( const point3D &p ) const
|
|
{
|
|
using std::max;
|
|
using std::abs;
|
|
return max( max( abs( p.x ) , abs( p.y ) ) , abs( p.z ) );
|
|
}
|
|
};
|
|
} } }
|
|
//]
|
|
|
|
std::ostream& operator<<( std::ostream &out , const point3D &p )
|
|
{
|
|
out << p.x << " " << p.y << " " << p.z;
|
|
return out;
|
|
}
|
|
|
|
//[point3D_main
|
|
const double sigma = 10.0;
|
|
const double R = 28.0;
|
|
const double b = 8.0 / 3.0;
|
|
|
|
void lorenz( const point3D &x , point3D &dxdt , const double t )
|
|
{
|
|
dxdt.x = sigma * ( x.y - x.x );
|
|
dxdt.y = R * x.x - x.y - x.x * x.z;
|
|
dxdt.z = -b * x.z + x.x * x.y;
|
|
}
|
|
|
|
using namespace boost::numeric::odeint;
|
|
|
|
int main()
|
|
{
|
|
|
|
point3D x( 10.0 , 5.0 , 5.0 );
|
|
// point type defines it's own operators -> use vector_space_algebra !
|
|
typedef runge_kutta_dopri5< point3D , double , point3D ,
|
|
double , vector_space_algebra > stepper;
|
|
int steps = integrate_adaptive( make_controlled<stepper>( 1E-10 , 1E-10 ) , lorenz , x ,
|
|
0.0 , 10.0 , 0.1 );
|
|
std::cout << x << std::endl;
|
|
std::cout << "steps: " << steps << std::endl;
|
|
}
|
|
//]
|