odeint - solving an ODE with time-dependent parameter numerically
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. 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.). This strategy is different and mathematically more accurate than assuming parameters are constant over a time step like a multistage method. Would anyone please help me to solve this problem using boost? Thank you in advance. With Best Regards, Arijit Hazra
On 1/19/2014, 5:55 PM, 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.
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.).
This strategy is different and mathematically more accurate than assuming parameters are constant over a time step like a multistage method.
Would anyone please help me to solve this problem using boost? Thank you in advance.
With Best Regards,
Arijit Hazra
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
Arijit, I'm not going to answer this from a Boost perspective (Boost mods jump in if need be), I'm going to answer this based on you asking this question over in SciComp at Stackexchange. I think you're looking for maths help, not necessarily Boost help. Having seen that other question and the answers, what are you looking for exactly? In your question, you said "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.". That's correct, you can. If a, b and c have widely-varying values (several orders of magnitude ratios between them) you have a stiff system. 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. 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? If you just need an answer you could use Maple, for example, and not have to program anything. You can even do it in Excel. Damien
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,
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; // ...
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
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
Hth, 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. 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@googlemail.com> 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@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
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.
Yes, then the RHS function (ode in the upper case) is passed to the stepper or integrate functions. Just take a look in the docs if this is not clear.^
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.
Do you need to exchange the functions for your system? In my suggestion I just created a function for coefficient a and and call this function from the rhs function. Otherwise you can make template paramters of it: template< typename CoefA , typename CoefB > struct ode { CoefA m_coef_a; CoefB m_coef_b; ode( CoefA coef_a , CoefB coef_b ) : m_coef_a( coef_a ) , m_coef_b( coef_b ) {} void operator()( state_type const& x , state_type &dxdt , double t ) const { dxdt[0] = m_coef_a( t ) * y[1] ... ... } }; maybe even with a nice make_ode function that you do not need to explicitly pass the template parameter template< typename CoefA , typename CoefB > ode< CoefA , CoefB > make_ode( CoefA coef_a , CoefB coef_b ) { return ode< CoefA , CoefB >( coef_a , coef_b ); } you can use this factory like auto my_ode = make_ode( concrete_a , concrete_b );
participants (4)
-
arijit hazra -
Arijit Hazra
-
Damien Hocking -
Karsten Ahnert