Add your own alternative version
Stats
137.4K views 2.3K downloads 96 bookmarked
Posted
9 Nov 2009

Comments and Discussions



Thank you, this will definitely be useful somehow





still think the same thing I thought in October and nothing's changed






Excellent stuff  well worth a 5+





Hi I'm running some code and I'm trying to use this library to get the code to run quicker. Basically the problem is, when I use my own RK4 (no error estimation, adaptive step size etc, just the basic thing), the code runs much quicker than if I implement the following:
runge_kutta4< state_type > stepper;
while(SRparameters[0] < 1) {
stepper.do_step(ODE, SRparameters, N, dN);
}
void ODE(const state_type &x, state_type &dxdN, const double N){
dxdN[0] = F(x[0], x[1]);
for(int l = 1; l < L; l++)dxdN[l] = G(x[0], x[1], x[l], x[l+1], l);
dxdN[L] = G(x[0], x[1], x[L], 0, L);
}
The difference is by a factor of 2 or 3 although its hard to say as I have suppressed a lot of what goes on. I'm not too sure but could the difference caused by the program calling from many different library's?
I was hoping to implement one of the more sophisticated methods but I might be doing something wrong to hinder the program. I'm compiling under Visual Studio 2010 if that makes any difference.





Sorry just realised I posted in the old version :(





Hi, its difficult to say whats going on. Can you post your rk4 method? Stupid question: Do you build the code in release mode?





Hi there thankyou for the quick reply.
void RK4(vector<double> &SR, double dt){
K1[0] = F(SR[0], SR[1])*dt;
for(int l = 1; l < L; l++){
K1[l] = G(SR[0], SR[1], SR[l], SR[l+1],l)*dt;
}
K1[L] = G(SR[0], SR[1], SR[L], 0.0, L)*dt;
K2[0] = F(SR[0] + 0.5*K1[0], SR[1] + 0.5*K1[1])*dt;
for(int l = 1; l < L; l++){
K2[l] = G(SR[0] + 0.5*K1[0], SR[1] + 0.5*K1[1], SR[l] + 0.5*K1[l], SR[l+1] + 0.5*K1[l+1],l)*dt;
}
K2[L] = G(SR[0] + 0.5*K1[0], SR[1] + 0.5*K1[1], SR[L] + 0.5*K1[L], 0.0,L)*dt;
K3[0] = F(SR[0] + 0.5*K2[0], SR[1] + 0.5*K2[1])*dt;
for(int l = 1; l < L; l++){
K3[l] = G(SR[0] + 0.5*K2[0], SR[1] + 0.5*K2[1], SR[l] + 0.5*K2[l], SR[l+1] + 0.5*K2[l+1],l)*dt;
}
K3[L] = G(SR[0] + 0.5*K2[0], SR[1] + 0.5*K2[1], SR[L] + 0.5*K2[L], 0.0,L)*dt;
K4[0] = F(SR[0] + K3[0], SR[1] + K3[1])*dt;
for(int l = 1; l < L; l++){
K4[l] = G(SR[0] + K3[0], SR[1] + K3[1], SR[l] + K3[l], SR[l+1] + K3[l+1],l)*dt;
}
K4[L] = G(SR[0] + K3[0], SR[1] + K3[1], SR[L] + K3[L], 0.0,L)*dt;
for(int l = 0; l < L+1; l++){
SR[l] = SR[l] + (1.0)/(6.0)*(K1[l] + 2*K2[l] + 2*K3[l] + K4[l]);
}
}
So far I've been running it in debug mode.





Ok I changed it to release mode and now it runs about 200x faster than my RK4. Thanks for the tip!





Ok, thats fine. I think your code should also be faster when compiled in release mode.





Hi sorry to bother you again but I'm having some problems when trying to get the code to integrate backwards in time and redefining time isn't really an option.
Basically I'm calling
integrate_adaptive(controlled_stepper, ODEpot, Observables, NvaluesSR[i], NvaluesSR[i1], dN, push_back_state_and_time(ObservableValues, NvaluesObs));
where NvaluesSR[i] > NvaluesSR[i1] and dN = 0.01. NvaluesSR are the time values push_back_state_and_time gave out earlier in the code when integrating (some different ODE's) forward in time.
Anyway at first the code wasn't doing anything, updating N or calling ODEpot. So I had a look in integrate_adaptive.hpp and figured the problem was caused by the loop (line 87)
while(start_time < end_time)
{
...
if( ( start_time + dt ) > end_time ){...}
...
}
which obviously wasn't running in my code as start_time > end_time. So to try and get it to work I just copied/pasted the same loop straight after but with < and > swapped. Anyway the problem seems to be that when a negative dt is passed through m_max_rel_error (controlled_runge_kutta.hpp*, line 238) becomes infinite and either the program grinds to a halt as it reduces dt to ridiculously tiny values or it crashes. I tried figuring out how to fix it but I ran into a deadend. Is this an easy fix? Any help would be greatly appreciated.





Hi, at the moment odeint is not able to work with negative time steps. But, if your system is autonomous, hence it does not explicitly depend on t, you can simply multiply the derivative with 1 within the evaluation of your ODE and it will integrate backward in time. Hope this helps.





Hi there, unfortunately the system does depend on t and the t dependence comes from integrating a function earlier in the code. So I only have the t dependance at discrete points (interpolating between them). I did try that before but, naturally it didn't work. Is there anything else I can?
Thankyou for the reply.





You can also introduce the time as a dependent variable with the trivial dynamic dt/dt = 1, or for the backward integration dt/dt = 1. This enlarges your system by one dimension, and you also might loose a bit of performance but this should be not so much. For example, the timedependent harmonic oscillator changes from
typedef array< double , 2 > state_type;
void harm_osc( const state_type &x , state_type &dxdt , double t )
{
dxdt[0] = x[1];
dxdt[1] =  x[0] + cos( omega * t );
}
to
typedef array< double , 3 > state_type;
void harm_osc( const state_type &x , state_type &dxdt , double t )
{
dxdt[0] = x[1];
dxdt[1] =  x[0] + cos( omega * x[2] );
dxdt[2] = 1.0;
}
You can also introduce the direction by one flag:
typedef array< double , 3 > state_type;
struct harm_osc
{
bool direction;
harm_osc( bool direction_ ) : direction( direction_ ) { }
void operator()( const state_type &x , state_type &dxdt , double t ) const
{
dxdt[0] = ( direction ) ? x[1] : x[1] ;
dxdt[1] = ( direction ) ?  x[0] + cos( omega * x[2] ) : x[0]  cos( omega * x[2] );
dxdt[2] = ( direction ) ? 1.0 : 1.0 ;
}
};
Hope this helps.






And your better suggestion is?





"Solving ordinary differential equations in C++ using BOOST on Linux"





That is definitely wrong, you can compile odeint on windows.





Thank you for this very much needed library.
Can you please provide example code for an adaptive integrator?
I try to use:
integrate_adaptive( stepper , lorenzo , x , 0.0 , 10.0 , dt );
but I get a compiler error that in essence says:
stepper_rk5_ck has no member named ‘try_step’
Just a very quick line or two hint on how to do this will be very helpful.
Thanks.





Ok, here is an example:
typedef stepper_rk5_ck< state_type > stepper_type;
controlled_stepper_standard< stepper_type >
rk5_controlled( 1.0e6 , 1.0e7 , 1.0 , 1.0 );
size_t steps = integrate_adaptive( rk5_controlled , lorenz , x , 0.0 , 10.0 , 1.0e4 ,
cout << _1 << tab << _2[0] << tab << _2[1] << "\n" )
For adaptive stepsize integration you need a controlled stepper which tries to perform one step. This controlled stepper typically uses an error stepper as shown in the above expample. You can also take a look in the examples of the development branch sandbox/odeint or the current documentation.
I will update the article very soon and include a section about adaptive stepsize integration.
Edit: Small fix in the example.
modified on Friday, May 14, 2010 2:57 AM





Thank you. The code snippet you provided works.
However, when I change
typedef stepper_rk5_ck< state_type > stepper_type;
to
typedef stepper_rk78_fehlberg< state_type > stepper_type;
I get the longwinded error pasted below that basically says that rk78_fehlberg does not have function do_step() .
Is this a bug with my code or is controlled stepping not implemented in rk78_fehlberg yet?
Thanks.
/Users/philipadmin/Packages/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp: In member function 'boost::numeric::odeint::controlled_step_result boost::numeric::odeint::controlled_stepper_standard<ErrorStepper>::try_step(DynamicalSystem&, typename ErrorStepper::container_type&, const typename ErrorStepper::container_type&, typename ErrorStepper::time_type&, typename ErrorStepper::time_type&) [with DynamicalSystem = VehicleDynamics, ErrorStepper = boost::numeric::odeint::stepper_rk78_fehlberg<std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::container_traits<std::vector<double, std::allocator<double> > > >]':
/Users/philipadmin/Packages/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp:199: instantiated from 'boost::numeric::odeint::controlled_step_result boost::numeric::odeint::controlled_stepper_standard<ErrorStepper>::try_step(DynamicalSystem&, typename ErrorStepper::container_type&, typename ErrorStepper::time_type&, typename ErrorStepper::time_type&) [with DynamicalSystem = VehicleDynamics, ErrorStepper = boost::numeric::odeint::stepper_rk78_fehlberg<std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::container_traits<std::vector<double, std::allocator<double> > > >]'
/Users/philipadmin/Packages/odeint/boost/numeric/odeint/integrator_adaptive_stepsize.hpp:54: instantiated from 'size_t boost::numeric::odeint::integrate_adaptive(ControlledStepper&, DynamicalSystem&, typename ControlledStepper::container_type&, typename ControlledStepper::time_type, typename ControlledStepper::time_type, typename ControlledStepper::time_type, Observer&) [with ControlledStepper = boost::numeric::odeint::controlled_stepper_standard<boost::numeric::odeint::stepper_rk78_fehlberg<std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::container_traits<std::vector<double, std::allocator<double> > > > >, DynamicalSystem = VehicleDynamics, Observer = void ()(double, std::vector<double, std::allocator<double> >&, VehicleDynamics&)]'
/Users/philipadmin/Packages/odeint/boost/numeric/odeint/integrator_adaptive_stepsize.hpp:87: instantiated from 'size_t boost::numeric::odeint::integrate_adaptive(ControlledStepper&, DynamicalSystem&, typename ControlledStepper::container_type&, typename ControlledStepper::time_type, typename ControlledStepper::time_type, typename ControlledStepper::time_type) [with ControlledStepper = boost::numeric::odeint::controlled_stepper_standard<boost::numeric::odeint::stepper_rk78_fehlberg<std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::container_traits<std::vector<double, std::allocator<double> > > > >, DynamicalSystem = VehicleDynamics]'
../src/TrueVehicle.cpp:34: instantiated from here
/Users/philipadmin/Packages/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp:158: error: no matching function for call to 'boost::numeric::odeint::stepper_rk78_fehlberg<std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::container_traits<std::vector<double, std::allocator<double> > > >::do_step(VehicleDynamics&, std::vector<double, std::allocator<double> >&, const std::vector<double, std::allocator<double> >&, double&, double&, std::vector<double, std::allocator<double> >&)'
make: *** [src/TrueVehicle.o] Error 1





Sry, this is bug. Thank you for finding them.
The bug can be easily fixed, if you change line 365 in stepper_rk78_fehlberg.hpp into const container_type &dxdt , . A simple add of the const specifier fixes this problem. The full definition of the function should look like:
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
container_type &x ,
const container_type &dxdt , time_type t ,
time_type dt ,
container_type &xerr )
{
...
}
Let me know if this works.





Thanks for looking into this.
I inserted const as you described and the code now compiles and runs.
I am concerned though that the integration results are incorrect.
When I use integrate_const I get the same results as MATLAB ode45 integrator.
However, when I use integrate_adaptive I get integration results that very clearly not correct, even when I crank down the tolerances.
I was also surprised to see that the number of steps taken is very small (six or seven) when I would expect at least 100 to get the accuracy I was asking for.
Would you mind doing a quick check on your side to ensure that adaptive integration is providing you accurate results.
Thanks.





Ok, I have checked it with Lorenz system and everything looks fine. What exactly have you done?
It is a bit difficult to check the results. We could use a system where we know the exact solution like the harmonic oscillator and then compare the results.
Edit: I have integrated the harmonic oscillator with adaptive stepsize integration and everything fits.
modified on Friday, May 14, 2010 3:50 AM





Adaptive step sizing appears to give spurious results. Hopefully, you can help me correct this.
Below, I have pasted code that compares constant odeint, adaptive odeint, and MATLAB ode45 (adaptive RK45) for the circular restricted threebody problem.
The constant step integrator and Matlab solution converge, but adaptive integration is out in left field somewhere and does not appear to converge as we tighten down error tolerances.
The output from the code follows:
Using constant step integrator with step size 1e06 : (May take 10 to 20 seconds ...)
0.133066 1.11213 0.134143 1.93131 0.383882 0.118057
Using adaptive step integrator (tolerances: 1e12 , 1e13 ) :
Used 21 steps.
0.301733 1.064 0.123245 1.84121 0.692931 0.126202
For reference, the MATLAB ode45 (tolerances: 1e6, 1e7) solution for this system is:
0.133066 1.11213 0.134143 1.93131 0.383882 0.118057
Here is the code to produce these results:
#include <iostream>
#include <iomanip>
#include <boost/numeric/odeint.hpp>
using namespace std;
using namespace boost::numeric::odeint;
typedef vector< double > state_type;
typedef stepper_rk78_fehlberg< state_type > stepper_type;
const double mu = 0.0123;
class DynamicalSystem {
public:
void operator()(state_type& q, state_type& dXdt, double t) {
double x = q[0];
double y = q[1];
double z = q[2];
double xp = q[3];
double yp = q[4];
double zp = q[5];
double invdE = 1.0/sqrt((x+mu)*(x+mu)+y*y+z*z);
double invdE3 = invdE*invdE*invdE;
double invdM = 1.0/sqrt((x1+mu)*(x1+mu)+y*y+z*z);
double invdM3 = invdM*invdM*invdM;
dXdt[0] = xp;
dXdt[1] = yp;
dXdt[2] = zp;
dXdt[3] = 2*yp + x  (1mu)*(x+mu)*invdE3  mu*(x1+mu)*invdM3;
dXdt[4] = 2*xp + y  (1mu)*y*invdE3  mu*y*invdM3;
dXdt[5] = (1.0mu)*invdE3*z  mu*invdM3*z;
}
};
int main( int argc , char **argv )
{
const double dt = 0.000001;
const double reltol = 1.0e12;
const double abstol = 1.0e13;
const int NX = 6;
state_type q0(NX,0.0);
q0[0] = 1.2;
q0[2] = 0.2;
q0[4] = 2.0;
double t0 = 0.0;
double tf = 1.0;
stepper_rk78_fehlberg< state_type > my_const_stepper;
my_const_stepper.adjust_size( q0 );
controlled_stepper_standard< stepper_type > my_controlled_stepper( reltol , abstol , 1.0 , 1.0 );
my_controlled_stepper.adjust_size( q0 );
DynamicalSystem my_dyn;
cout << "\nUsing constant step integrator with step size " << dt << " : (May take 10 to 20 seconds ...)" << endl;
state_type xc(q0);
integrate_const( my_const_stepper , my_dyn , xc , t0 ,tf , dt );
for(int ii=0;ii<xc.size();++ii) cout << setw(12) << xc[ii] << flush;
cout << endl;
cout << "\nUsing adaptive step integrator (tolerances: " << reltol << " , " << abstol << " ) :" << endl;
state_type xa(q0);
size_t steps = integrate_adaptive( my_controlled_stepper , my_dyn , xa , t0 , tf , dt );
cout << "Used " << steps << " steps." << endl;
for(int ii=0;ii<xa.size();++ii) cout << setw(12) << xa[ii] << flush;
cout << endl;
cout << "\nFor reference, the MATLAB ode45 (tolerances: 1e6, 1e7) solution for this system is:" << endl;
cout << setw(12) << 0.133066 << setw(12) << 1.112127 << setw(12) << 0.134143 << setw(12) << 1.931314 << setw(12) << 0.383882 << setw(12) << 0.118057 << endl << endl;
return 0;
}





Ok, I figured out what happens. The adaptive integration routine does not stop exactly at t = tf . It performs the stepping as long as t < tf . If the time of one step is larger then tf it stops. So, the result you see has simply another timestep, larger then tf .
Of course, this is an unexpected behavior and has to be fixed. It should stop exactly at t = tend . I will do this very soon.





Hi headmyshoulder,
[bug fix]
goto "integrator_adaptive_stepsize.hpp" line 53
change
if( (end_time  t) < dt ) dt = end_time  t;
to
if( (end_time  t) < dt_ ) dt_ = end_time  t;
dt_ changes inside the loop, dt doesn't.
Thx for the codsnip safariman, it really did save me a lot of time, as well as odeint.
New output for the example from safariman
Using constant step integrator with step size 1e06 : (May take 10 to 20 seconds ...)
0.133066 1.11213 0.134143 1.93131 0.383882 0.118057
Using adaptive step integrator (tolerances: 1e12 , 1e13 ) :
Used 21 steps.
0.133066 1.11213 0.134143 1.93131 0.383882 0.118057
For reference, the MATLAB ode45 (tolerances: 1e6, 1e7) solution for this system is:
0.133066 1.11213 0.134143 1.93131 0.383882 0.118057
cheers,
Gregor.






headmyshoulder wrote: rk5_controlled( 10e6 , 1.0e7 , 1.0 , 1.0 );
should be
headmyshoulder wrote: rk5_controlled( 1.0e6 , 1.0e7 , 1.0 , 1.0 );





Thanks, I have updated the example.





Is there any way to compile this using Visual C++? I am not that profficient in Linux..please, it would be of great help





Hi, there is a comment on how to compile odeint under VS2008: [1]. Nevertheless, I will add a compile section for Visual Studio very soon.





Is there some way that a wrapper could be placed around this for use in C#. This may be a really dense question but it's been a long long time since I used C++ at all.
Thanks
Ken





I don't have any experience with C#, so I can't answer your question. But I have the feeling, that it is possible to call odeint from C#. A problem might be, that odeint uses heavily templates and as far as I know C# does not support them. Maybe one has to specialize them with a general container like the STL vectors.






I think doing the exactly the same in C# is difficult. But, if you "only" want to solve an ODE you can implement the Runge Kutta method. In C++ the minimal version looks like
void deriv( const double *x , double *dxdt , double t )
{
}
void rk4( double *x , size_t n , double t , double dt )
{
double dh = 0.5 * dt , d6= dt / 6.0;
double th = t + dh;
double dxdt[n] , dxm[n] , dxt[n] , xt[n];
deriv( x , dxdt , t );
for( size_t i=0 ; i<n ; i++ ) xt[i] = x[i] + dh * dxdt[i];
deriv( xt , dxt , th );
for( size_t i=0 ; i<n ; i++ ) xt[i] = x[i] + dh * dxt[i];
deriv( xt , dxm , th );
for( size_t i=0 ; i<n ; i++ )
{
xt[i] = x[i] + dt * dxm[i];
dxm[i] += dxt[i];
}
deriv( xt , dxt , t+dt );
for( size_t i=0 ; i<n ; i++ )
x[i] += d6 * ( dxdt[i] + dxt[i] + 2.0 * dxm[i] );
}
size_t n = 3; double x[n];
const double dt = 0.01;
double t = 0.0;
for( size_t ii=0 ; ii<10000 ; ++ii , t+=dt )
{
rk4( x , 3 , t , dt ); }
I think, you can easily integrate the above lines of code in C#. Hope, this helps.
modified on Monday, March 29, 2010 3:37 PM





There is a typo in the namespace in the text. It should read: boost::numeric::odeint. It is shown correctly in the code.
By the way, if anyone is interested in compiling under Ms VisualStudio, go to the main Boost pages and find section Getting started. It has a nice stepbystep HowTo build Boost on Windows system.
I tried in Visual Studio 2008 and needed to do following (downloaded full Boost depository):
 Add additional Include directory > path to the Boost_Root_Dir
 Create a New Project from existing code in File menu
 chosed a C++ project
 Project file directory is where the odeint.zip was unwrapped
 Build the project using "Visual Studio"
 set type of the project to "Console application"
 ignore the rest and click finish
and voiala... click Run.
btw. The first time I had to do following modification to avoid an error message:
in file "container_traits_tr1_array.hpp"
find the line: #include <tr1/array>
and change it to: #include <boost/tr1/array.hpp>
modified on Wednesday, March 24, 2010 12:33 PM





Thanks you very much for that hint and the MSVC guide. Unfortunately, MSVC and gcc have different locations for the tr1 headers, therefore you had to uninclude the "container_traits_tr1_array.hpp" line. I do not know, how one can include these headers consistently. Maybe anyone else knows?





Hi. First of all, thanks for such an elegant and apparently simple method to solve ODEs numerically. Only one question: Could you add a little bit about how to compile it? I am not very proficient in C++, and I am getting many errors trying to compile it with MingW. A simple paragraph with 1. do this, 2. do that, 3... etc to compile one of your examples could help very much. In essence, many of us just need a model to imitate. Thanks!





Ok, thanks for this hint. I will do this very soon. The whole library is headeronly and you only have to include the odeint.hpp. Nevertheless, a Makefile is provided with the sources, doing the job of compiling the examples, at least at most unixoid os.
Alternativly you can compile everything by yourself and you only have to provide the current working directory as a header search path. With g++ or gcc the command is the following:
g++ I. expample.cc
For lorenz_integrate_constant_step.cpp you also need the boost library. Here, the command is
g++ I. I$(BOOST_ROOT) expample.cc





Thanks so much for your response! I had tried the first option before. This is what I get when trying to compile from a folder in which I unzipped all files you provided:
C:\odeint>g++ I. lorenz_simple.cpp
lorenz_simple.cpp:20:21: tr1/array: No such file or directory
In file included from ./odeint/euler.hpp:17,
from ./odeint.hpp:12,
from lorenz_simple.cpp:22:
./odeint/resizer.hpp:47: error: `std::tr1' has not been declared
./odeint/resizer.hpp:47: error: `array' was not declared in this scope
./odeint/resizer.hpp:47: error: wrong number of template arguments (2, should be
1)
./odeint/resizer.hpp:18: error: provided for `template<class Container> class od
eint::resizer'
./odeint/resizer.hpp:47: error: expected unqualifiedid before '>' token
lorenz_simple.cpp:30: error: `std::tr1' has not been declared
lorenz_simple.cpp:30: error: expected initializer before '<' token
lorenz_simple.cpp:36: error: variable or field `lorenz' declared void
lorenz_simple.cpp:36: error: `state_type' was not declared in this scope
lorenz_simple.cpp:36: error: `x' was not declared in this scope
lorenz_simple.cpp:36: error: `state_type' was not declared in this scope
lorenz_simple.cpp:36: error: `dxdt' was not declared in this scope
lorenz_simple.cpp:36: error: expected primaryexpression before "double"
lorenz_simple.cpp:37: error: initializer expression list treated as compound exp
ression
lorenz_simple.cpp:37: error: expected `,' or `;' before '{' token
lorenz_simple.cpp:46: error: `state_type' has not been declared
lorenz_simple.cpp:46: error: `state_type' has not been declared
lorenz_simple.cpp:47: error: ISO C++ forbids declaration of `x' with no type
lorenz_simple.cpp:47: error: ISO C++ forbids declaration of `dxdt' with no type
lorenz_simple.cpp: In member function `void lorenz_class::operator()(int&, int&,
double)':
lorenz_simple.cpp:48: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp:48: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp:48: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp:49: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp:49: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp:49: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp:49: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp:49: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp:50: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp:50: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp:50: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp:50: error: invalid types `int[int]' for array subscript
lorenz_simple.cpp: In function `int main(int, char**)':
lorenz_simple.cpp:58: error: `state_type' was not declared in this scope
lorenz_simple.cpp:58: error: expected `;' before "x"
lorenz_simple.cpp:59: error: `state_type' cannot appear in a constantexpression
lorenz_simple.cpp:59: error: template argument 1 is invalid
lorenz_simple.cpp:59: error: template argument 3 is invalid
lorenz_simple.cpp:59: error: invalid type in declaration before ';' token
lorenz_simple.cpp:64: error: request for member `next_step' in `stepper', which
is of nonclass type `int'
lorenz_simple.cpp:64: error: `x' was not declared in this scope
Any ideas? Thanks in advance! BTW, running Win7, MingW 5.1.6





Ok, your platform does not support tr1 headers [1], or at least does not find them. I have no adhoc solution for this problem right now, but will change this soon. In the meantime you can do is the following:
1. remove the resizer< array > specialization in odeint/resizer.hpp
2. change typedef array< double , n > state_type to typedef vector< double > state_type
3. allocate enough memory for the vectors.
One examples shows, that arrays are faster due to their fixed size.





Thanks for that hint. Now there seems to be a smaller problem, which I am not sure how to solve:
C:\odeint>g++ I. lorenz_simple.cpp
lorenz_simple.cpp: In function `int main(int, char**)':
lorenz_simple.cpp:59: error: `x' must be initialized by constructor, not by `{...}'
Any additional hint? I hope I am no offending you with my lack of expertise in C++. The bright side is that finding and solving these difficulties will make your library more accessible to people like me: competent mathematicians but mediocre programmers.
Thanks!





This is still an relict from tr1/array . If you have changed state_type to vector (what I suppose, since no other errors occur) you only have to replace the line
state_type x = {{ 1.0 , 0.0 , 0.0 }}; by
state_type x(3);
x[0] = 1.0;
x[1] = 0.0;
x[2] = 0.0
This initializes a vector of length 3 with the values (1,0,0). The double braces can only be used to initialize tr1/array 's. Hope this helps.





It worked!!! THANKS for your help. I will solve my problem, document a stepbystep solution in my web page, and post the link here relatively soon. Hopefully that will help others too. I had to quit on Matlab because the time it took to produce a bifurcation manifold was unacceptable. I was looking simultaneously at the classic DLSODE and a couple other options. I could produce results first with with this library. GREAT JOB! It is soooo simple to manipulate the dynamical system I am studying... Please keep improving this library and expanding it. I have plenty of research models for testing and will be glad to contribute (my focus is on analysis, but numerics is unavoidable).
All best,
J.





When I took Numerical Analysis which included error analysis of various techniques for numerical solution of ODES, From Atkinson "An Introduction to Numerical Analysis". The RungeKutta Fehlberg 4 order combined with the extra evaluation of the 5th order gave an automatic error analysis and thus allowed for automatic adjustment of the step size used in the solution of the ODE's. This method could be used for many many ODE's without fear of loss of error, since the epsilon for the computer 1.0+Epsilon > 1.0 where Epsilon is as small as possible.
The RKF45 method was also included in the "Computer Methods for Mathematical Computations" by Forsythe Malcom and Moler.
My point is that you should never be using a straight RungeKutta method, but should always be using the Fehlberg variation to improve error analysis and allow your program to adjust the step size, according to the error results by comparing the 4th and 5th order methods. Of course if the objective function is prohibitively expensive, then this method may be expensive, in which case other methods such as collocation, or spectral methods may be better. For a general ODE solver the RKF45 method with step adjustment simply has stood the test of time.





Error analysis and stepsize control is already included in the current version of odeint. Up to now, the RK45 method is missing, but a method which estimates the error by calculating two steps with dt/2 is present. RK45 and RK78 will be included very soon. Nevertheless, in many cases one does not need adaptive stepsize control, think for example on creating artificial time series from an ODE.





I can assure you that we are aware of the standard numerical methods to solve ode's. The article is focusing more on the structural part of the code and the basic interfaces to some (abstract) step and error adaption routines. We are currently implementing more sophisticated integration methods, like the rk45 with automatic error evaluation (I've just commited the code for that). RK89 will follow as well as BulirschStoer methods, I think.







General News Suggestion Question Bug Answer Joke Praise Rant Admin Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

