123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169 |
- /*
- * fpu.cpp
- *
- * This example demonstrates how one can use odeint to solve the Fermi-Pasta-Ulam system.
- * Created on: July 13, 2011
- *
- * Copyright 2011-2012 Karsten Ahnert
- * Copyright 2011 Mario Mulansky
- * 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)
- */
- #include <iostream>
- #include <numeric>
- #include <cmath>
- #include <vector>
- #include <boost/numeric/odeint.hpp>
- #ifndef M_PI //not there on windows
- #define M_PI 3.1415927 //...
- #endif
- using namespace std;
- using namespace boost::numeric::odeint;
- //[ fpu_system_function
- typedef vector< double > container_type;
- struct fpu
- {
- const double m_beta;
- fpu( const double beta = 1.0 ) : m_beta( beta ) { }
- // system function defining the ODE
- void operator()( const container_type &q , container_type &dpdt ) const
- {
- size_t n = q.size();
- double tmp = q[0] - 0.0;
- double tmp2 = tmp + m_beta * tmp * tmp * tmp;
- dpdt[0] = -tmp2;
- for( size_t i=0 ; i<n-1 ; ++i )
- {
- tmp = q[i+1] - q[i];
- tmp2 = tmp + m_beta * tmp * tmp * tmp;
- dpdt[i] += tmp2;
- dpdt[i+1] = -tmp2;
- }
- tmp = - q[n-1];
- tmp2 = tmp + m_beta * tmp * tmp * tmp;
- dpdt[n-1] += tmp2;
- }
- // calculates the energy of the system
- double energy( const container_type &q , const container_type &p ) const
- {
- // ...
- //<-
- double energy = 0.0;
- size_t n = q.size();
- double tmp = q[0];
- energy += 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
- for( size_t i=0 ; i<n-1 ; ++i )
- {
- tmp = q[i+1] - q[i];
- energy += 0.5 * ( p[i] * p[i] + tmp * tmp ) + 0.25 * m_beta * tmp * tmp * tmp * tmp;
- }
- energy += 0.5 * p[n-1] * p[n-1];
- tmp = q[n-1];
- energy += 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
- return energy;
- //->
- }
- // calculates the local energy of the system
- void local_energy( const container_type &q , const container_type &p , container_type &e ) const
- {
- // ...
- //<-
- size_t n = q.size();
- double tmp = q[0];
- double tmp2 = 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
- e[0] = tmp2;
- for( size_t i=0 ; i<n-1 ; ++i )
- {
- tmp = q[i+1] - q[i];
- tmp2 = 0.25 * tmp * tmp + 0.125 * m_beta * tmp * tmp * tmp * tmp;
- e[i] += 0.5 * p[i] * p[i] + tmp2 ;
- e[i+1] = tmp2;
- }
- tmp = q[n-1];
- tmp2 = 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
- e[n-1] += 0.5 * p[n-1] * p[n-1] + tmp2;
- //->
- }
- };
- //]
- //[ fpu_observer
- struct streaming_observer
- {
- std::ostream& m_out;
- const fpu &m_fpu;
- size_t m_write_every;
- size_t m_count;
- streaming_observer( std::ostream &out , const fpu &f , size_t write_every = 100 )
- : m_out( out ) , m_fpu( f ) , m_write_every( write_every ) , m_count( 0 ) { }
- template< class State >
- void operator()( const State &x , double t )
- {
- if( ( m_count % m_write_every ) == 0 )
- {
- container_type &q = x.first;
- container_type &p = x.second;
- container_type energy( q.size() );
- m_fpu.local_energy( q , p , energy );
- for( size_t i=0 ; i<q.size() ; ++i )
- {
- m_out << t << "\t" << i << "\t" << q[i] << "\t" << p[i] << "\t" << energy[i] << "\n";
- }
- m_out << "\n";
- clog << t << "\t" << accumulate( energy.begin() , energy.end() , 0.0 ) << "\n";
- }
- ++m_count;
- }
- };
- //]
- int main( int argc , char **argv )
- {
- //[ fpu_integration
- const size_t n = 64;
- container_type q( n , 0.0 ) , p( n , 0.0 );
- for( size_t i=0 ; i<n ; ++i )
- {
- p[i] = 0.0;
- q[i] = 32.0 * sin( double( i + 1 ) / double( n + 1 ) * M_PI );
- }
- const double dt = 0.1;
- typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
- fpu fpu_instance( 8.0 );
- integrate_const( stepper_type() , fpu_instance ,
- make_pair( boost::ref( q ) , boost::ref( p ) ) ,
- 0.0 , 1000.0 , dt , streaming_observer( cout , fpu_instance , 10 ) );
- //]
- return 0;
- }
|