Re: [Boost-users] odeint - solving an ODE with time-dependent parameter numerically

Subject: Re: [Boost-users] odeint - solving an ODE with time-dependent parameter numerically
From: Karsten Ahnert (karsten.ahnert_at_[hidden])
Date: 2014-01-20 01:30:03

Hi,

On 01/20/2014 01:55 AM, Arijit Hazra wrote:
> Hello All,
> 
> I am trying to solve an ode with time-dependent parameters by using
> boost. I have just started using boost and I think boost library's
> odeint library is the best option for my problem.
> 
> Would anyone please explain me how I can solve a Initial Value Problem
> --a system of ODE --with a time-dependent parameter. e.g.
> 
> |   dy1/dt = a1(t) * y2 -b1(t)* y3 -c1(t) y1;
>    dy2/dt = a2(t) * y1 -b2(t)* y1 -c2(t) y3;
>    dy3/dt = a3(t) * y2 -b3(t)* y3 -c3(t) y1;
> |
> 
> y(0)=c; My idea so far after  going through the manuals and help files.
> is that for fixed time-step if I form a vector of a,b,c and assume a,b,c
> to be constant over that time-period I can solve using boost library.

How do a,b,c look like? Are it functions you can evaluate, like sin(),
cos(), .... or combinations of it? If this the case you can just write
down the equations like

double a1( double t ) { // return a1(t); }
double a2( double t ) { // return a1(t); }
// ...

void ode( state_type const &y , state_type &dydt , double t )
{
    dydt[0] = a1(t) * y[1] - b1(t) * y[2] - c1(t) * y[0];
    // ...
}

The steppers in odeint do not assume that a,b,c are constant over the
time step. If you use for example a multi-stage method like Runge-Kutta
it will evaluate a,b,c at different times during one time step. It also
works for controlled stepper with step size control. Ff course, the
stepper do not know what happens between the evaluations - this
inaccuracy is the price you have to pay when solving ODEs numerically,
but it should not be too large and can be controlled when using step
size adaption.

> But this is mathematically less accurate and computationally expensive 
> as far my understanding .For solving an ODE numerically
> /y/Ë™(/t/)=/f/(/y/(/t/),/t/) generally the problem is all about function
> evaluation for explicit methods and Jacobian evaluation for Implicit
> methods at specific points in /t/ and /y/.
> 
> e.g a step of explicit Euler is /yn/+1=/yn/+(/tn/+1−/tn/)⋅/f/(/yn/,/tn/).
> 
> For time-dependent parameter, treating it as part of /f/, the function
> evaluation (respectively, for implicit methods, also treat it as part of
> the Jacobian evaluation) is the standard procedure. A similar strategy
> applies to more complicated methods for solving ODEs (multistage methods
> such as Runge-Kutta, implicit methods for stiff systems, etc.).

Nearly all steppers pass the time explicitly to the ODE, see above. You
can safely use this parameter to solve a non-autonomous ODE. Of course,
this is exactly equivalent to solve the enhanced ODE with an additional
dimension for the time, in your example x1 = y1 , x2= y2 , x3 = y3 , x4 = t;

therefore

dx1/dt = a1(x4) * x2 - b1(x4) * x3 - c1(x4) * x1;
// ...

Hth,