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
fpufpu36 fpu( const double beta = 1.0 ) : m_beta( beta ) { }
37
38 // system function defining the ODE
operator ()fpu39 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
energyfpu58 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
local_energyfpu81 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
streaming_observerstreaming_observer114 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 >
operator ()streaming_observer118 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
main(int argc,char ** argv)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 }
170