Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Vector Runge-Kutta ODE solver with compilation?

• To: mathgroup at smc.vnet.net
• Subject: [mg116611] Re: Vector Runge-Kutta ODE solver with compilation?
• From: DmitryG <einschlag at gmail.com>
• Date: Mon, 21 Feb 2011 19:29:12 -0500 (EST)
• References: <201102210919.EAA20606@smc.vnet.net> <ijth45$mcp$1@smc.vnet.net>

On 21 Feb., 06:08, Oliver Ruebenkoenig <ruebe... at wolfram.com> wrote:
> DmitryG,
>
> one thing I forgot to mention is
>
> tutorial/NDSolvePlugIns
>
> that will show you how you can then use your own method and plug that into
> NDSolve.
>
> Oliver
>
>
>
>
>
>
>
> On Mon, 21 Feb 2011, DmitryG wrote:
> > Dear All,
>
> > Inspired by the new capability of Mathematica 8 to compile pieces of
> > user code with its own or external C compiler, I am trying to
> > implement the simple 4th-order Runge-Kutta solver with fixed step in a
> > vector form, that is, suitable for solving an arbitrary number of
> > equations.
>
> > The solver without compilation, the idea of which is taken from
> >http://www.jasoncantarella.com, is working fine:
>
> > (* Definition *)
> > RK4[f_,x0_,t0_,tMax_,n_] := Module[{h,K1,K2, K3,
> > K4,x=x0,SolList={{t0,x0}}},
>
> > h = (tMax - t0)/n;
>
> > Do[
> > t=t0+k h;
> > K1 = h f[t,x];
> > K2 = h f[t+(1/2) h, x +(1/2)  K1];
> > K3 = h f[t+(1/2)  h, x + (1/2)  K2];
> > K4 = h f[t+h, x + K3];
> > x=x + (1/6) K1 +(1/3) K2 +(1/3) K3 +(1/6) K4;
>
> > SolList=Append[SolList,{t,x}]
>
> > ,{k,1,n}];
> > SolList
> > ];
>
> > (* Testing for a system of two equations *)
> > F[t_,x_] := {x[[2]],1 - x[[1]] + 4/(x[[1]]^3)};
> > t0=AbsoluteTime[];
> > Solution=RK4[ F,{1.0,1.0},0.0,200.0,5000];
> > AbsoluteTime[]-t0
> > ListLinePlot[Take[Solution[[All,2]],100],PlotMarkers-
> >> Automatic,PlotStyle->{Red}]
>
> > Out[55]= 0.5468750
>
> > The output above is execution time. The original code uses NestList
> > instead of my Do cycle but the execution time is the same.
>
> > Unfortunately, the compiled version of this code does not work as
> > expected:
>
> > (* Definition *)
> > RK4Comp =
> >  Compile[{{f, _Real, 1}, {x0, _Real, 1}, {t0, _Real}, {tMax, _Real},
> > {n, _Integer}},
> >   Module[{h, K1, K2, K3, K4, SolList = {{t0, x0}}, x = x0, t},
>
> >    h = (tMax - t0)/n;
>
> >    Do[
> >     t = t0 + k h;
> >     K1 = h f[t, x];
> >     K2 = h f[t + (1/2) h, x + (1/2)  K1];
> >     K3 = h f[t + (1/2)  h, x + (1/2)  K2];
> >     K4 = h f[t + h, x + K3];
> >     x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
>
> >     SolList = Append[SolList, {t, x}]
> >     , {k, 1, n}];
> >    SolList
> >    ](*,CompilationTarget->"C"*)
> >   ]
>
> > (* Testing for a system of two ODEs *)
> > F[t_, x_] := {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)};
> > t0 = AbsoluteTime[];
> > Solution = RK4Comp[ F, {1.0, 1.0}, 0.0, 200.0, 5000];
> > AbsoluteTime[] - t0
> > ListLinePlot[Take[Solution[[All, 2]], 100], PlotMarkers -> Automatic,
> > PlotStyle -> {Red}]
>
> > During evaluation of In[57]:= CompiledFunction::cfta: Argument F at
> > position 1 should be a rank 1 tensor of machine-size real numbers. >>
>
> > Out[61]= 0.5312500
>
> > Mathematica complains and seemingly proceeds without compilation
> > because execution time is the same.
>
> > Anybody has an idea of what is happening and how it can be corected? I
> > believe the problem should be relevant for many.
>
> > Thank you for your time,
>
> > Dmitry

Oliver,

This looks nice and works but is very slow. What is the step size h in
the fixed-step RK4 algorithm? Here is my working version for the same
equations:

**************************************************************************

CRK4[]["Step"[rhs_, t_, h_, y_, yp_]] := Module[{k0, k1, k2, k3},
  k0 = h yp;
  k1 = h rhs[t + h/2, y + k0/2];
  k2 = h rhs[t + h/2, y + k1/2];
  k3 = h rhs[t + h, y + k2];
  {h, (k0 + 2  k1 + 2  k2 + k3)/6}]
CRK4[___]["DifferenceOrder"] := 4
CRK4[___]["StepMode"] := Fixed   (* What is "Fixed"?? Mathematica does
not know this. What is the step h?? *)

Equations = {x1'[t] == x2[t], x2'[t] == 1 - x1[t] + 4/(x1[t]^3)};
IniConds = {x1[0] == 1, x2[0] == 1};
tMax = 200;

tt0 = AbsoluteTime[];
fixed = NDSolve[Join[Equations, IniConds], {x1, x2}, {t, 0, tMax},
Method -> CRK4, MaxSteps -> 1000000]
AbsoluteTime[] - tt0

x1t[t_] := x1[t] /. fixed[[1]];
x2t[t_] := x2[t] /. fixed[[1]];

ParametricPlot[{x1t[t], x2t[t]}, {t, 0, tMax}]

Out[61]=3.0000000

***************************************************************************=
*

Here, for tMax=500, the execution time is 3.00 and the plot shows that
the precision is insufficient because of a too large h. For larger
tMax it becomes dramatic.

Best,

Dmitry

• References:
  • Vector Runge-Kutta ODE solver with compilation?
    • From: DmitryG <einschlag@gmail.com>