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: arijit hazra (mailtohazra_at_[hidden])
Date: 2014-01-21 08:10:12

> 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

my parameters are basically a1(t,params) ,a2(t,params), a3(t,params) are
piecewise continuous functions dependent on t and some other fixed
parameters e.g. it is sinc function in some interval,
constant in some interval of t and linear in some interval of t etc.

>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];
   > // ...
>}

I am little bit confused now how to incorporate the RHS of the ode system-I
guess I need to write RHS of the ode separately.Your suggestion would be
helpful.
e.g.

double a1(double t, double c1, double c2)
{
//
return a1;
}

similarly for a2 and a3. i am confused how I can pass the a1,b1,c1 and so
on in the ODE function and solve it.
Suggestions would be helpful.

>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,

>Different algorithms are derived to solve different ways, some assume a, b
and c are constant over the whole step, and >some require that you update
a, b and c during the substeps of the calculation (e.g. multistage).
Typically, algorithms >assume you are going to do the "right thing" and
that whenever they ask you to evaluate dydt, you're evaluating that >vector
with the most recent information.  You don't *have* to do that in a
multistage algorithm, but as the SciComp answer >suggested, if you do
assume your parameters are constant, you're effectively wasting the
advantages of a multistage >method.

Now, I understood the mathematical parts I need though it is far from
rigorous knowledge.Thanks for your cooperation.

>Is this homework?  Do you have to do this in C++ (which is what Boost is
written in), or can you choose any language?  >Do you actually have to
program it?

No this is not a homework problem. I am trying to build a code with solving
the system of ODE I mentioned as the first building block.I think odeint
library of boost is best for my purpose.

Thank you for your co-operation.

Arijit

On Mon, Jan 20, 2014 at 7:30 AM, Karsten Ahnert <
karsten.ahnert_at_[hidden]> wrote:

> 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,
>
>
> _______________________________________________
> Boost-users mailing list
> Boost-users_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/boost-users
>