libs/numeric/odeint/examples/fpu.cpp

   1 /*
   2  * fpu.cpp
   3  *
   4  * This example demonstrates how one can use odeint to solve the Fermi-Pasta-Ulam system.
   5
   6  *  Created on: July 13, 2011
   7  *
   8  * Copyright 2011-2012 Karsten Ahnert
   9  * Copyright 2011 Mario Mulansky
  10  * Distributed under the Boost Software License, Version 1.0. (See
  11  * accompanying file LICENSE_1_0.txt or copy at
  12  * http://www.boost.org/LICENSE_1_0.txt)
  13  */
  14
  15 #include <iostream>
  16 #include <numeric>
  17 #include <cmath>
  18 #include <vector>
  19
  20 #include <boost/numeric/odeint.hpp>
  21
  22 #ifndef M_PI //not there on windows
  23 #define M_PI 3.1415927 //...
  24 #endif
  25
  26 using namespace std;
  27 using namespace boost::numeric::odeint;
  28
  29 //[ fpu_system_function
  30 typedef vector< double > container_type;
  31
  32 struct fpu
  33 {
  34     const double m_beta;
  35
  36     fpu( const double beta = 1.0 ) : m_beta( beta ) { }
  37
  38     // system function defining the ODE
  39     void operator()( const container_type &q , container_type &dpdt ) const
  40     {
  41         size_t n = q.size();
  42         double tmp = q[0] - 0.0;
  43         double tmp2 = tmp + m_beta * tmp * tmp * tmp;
  44         dpdt[0] = -tmp2;
  45         for( size_t i=0 ; i<n-1 ; ++i )
  46         {
  47             tmp = q[i+1] - q[i];
  48             tmp2 = tmp + m_beta * tmp * tmp * tmp;
  49             dpdt[i] += tmp2;
  50             dpdt[i+1] = -tmp2;
  51         }
  52         tmp = - q[n-1];
  53         tmp2 = tmp + m_beta * tmp * tmp * tmp;
  54         dpdt[n-1] += tmp2;
  55     }
  56
  57     // calculates the energy of the system
  58     double energy( const container_type &q , const container_type &p ) const
  59     {
  60         // ...
  61         //<-
  62         double energy = 0.0;
  63         size_t n = q.size();
  64
  65         double tmp = q[0];
  66         energy += 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
  67         for( size_t i=0 ; i<n-1 ; ++i )
  68         {
  69             tmp = q[i+1] - q[i];
  70             energy += 0.5 * ( p[i] * p[i] + tmp * tmp ) + 0.25 * m_beta * tmp * tmp * tmp * tmp;
  71         }
  72         energy += 0.5 * p[n-1] * p[n-1];
  73         tmp = q[n-1];
  74         energy += 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
  75
  76         return energy;
  77         //->
  78     }
  79
  80     // calculates the local energy of the system
  81     void local_energy( const container_type &q , const container_type &p , container_type &e ) const
  82     {
  83         // ...
  84         //<-
  85         size_t n = q.size();
  86         double tmp = q[0];
  87         double tmp2 = 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
  88         e[0] = tmp2;
  89         for( size_t i=0 ; i<n-1 ; ++i )
  90         {
  91             tmp = q[i+1] - q[i];
  92             tmp2 = 0.25 * tmp * tmp + 0.125 * m_beta * tmp * tmp * tmp * tmp;
  93             e[i] += 0.5 * p[i] * p[i] + tmp2 ;
  94             e[i+1] = tmp2;
  95         }
  96         tmp = q[n-1];
  97         tmp2 = 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
  98         e[n-1] += 0.5 * p[n-1] * p[n-1] + tmp2;
  99         //->
 100     }
 101 };
 102 //]
 103
 104
 105
 106 //[ fpu_observer
 107 struct streaming_observer
 108 {
 109     std::ostream& m_out;
 110     const fpu &m_fpu;
 111     size_t m_write_every;
 112     size_t m_count;
 113
 114     streaming_observer( std::ostream &out , const fpu &f , size_t write_every = 100 )
 115     : m_out( out ) , m_fpu( f ) , m_write_every( write_every ) , m_count( 0 ) { }
 116
 117     template< class State >
 118     void operator()( const State &x , double t )
 119     {
 120         if( ( m_count % m_write_every ) == 0 )
 121         {
 122             container_type &q = x.first;
 123             container_type &p = x.second;
 124             container_type energy( q.size() );
 125             m_fpu.local_energy( q , p , energy );
 126             for( size_t i=0 ; i<q.size() ; ++i )
 127             {
 128                 m_out << t << "\t" << i << "\t" << q[i] << "\t" << p[i] << "\t" << energy[i] << "\n";
 129             }
 130             m_out << "\n";
 131             clog << t << "\t" << accumulate( energy.begin() , energy.end() , 0.0 ) << "\n";
 132         }
 133         ++m_count;
 134     }
 135 };
 136 //]
 137
 138
 139
 140
 141
 142
 143
 144
 145 int main( int argc , char **argv )
 146 {
 147     //[ fpu_integration
 148     const size_t n = 64;
 149     container_type q( n , 0.0 ) , p( n , 0.0 );
 150
 151     for( size_t i=0 ; i<n ; ++i )
 152     {
 153         p[i] = 0.0;
 154         q[i] = 32.0 * sin( double( i + 1 ) / double( n + 1 ) * M_PI );
 155     }
 156
 157
 158     const double dt = 0.1;
 159
 160     typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
 161     fpu fpu_instance( 8.0 );
 162
 163     integrate_const( stepper_type() , fpu_instance ,
 164             make_pair( boost::ref( q ) , boost::ref( p ) ) ,
 165             0.0 , 1000.0 , dt , streaming_observer( cout , fpu_instance , 10 ) );
 166     //]
 167
 168     return 0;
 169 }