       Re: Vector Runge-Kutta ODE solver with compilation?

• To: mathgroup at smc.vnet.net
• Subject: [mg116650] Re: Vector Runge-Kutta ODE solver with compilation?
• From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
• Date: Tue, 22 Feb 2011 06:25:34 -0500 (EST)

```On Mon, 21 Feb 2011, DmitryG wrote:

> 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[],1 - x[] + 4/(x[]^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= 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[], 1 - x[] + 4/(x[]^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:= CompiledFunction::cfta: Argument F at
>>> position 1 should be a rank 1 tensor of machine-size real numbers. >>
>>
>>> Out= 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
>

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 == 1, x2 == 1};
> tMax = 200;
>
> tt0 = AbsoluteTime[];
> fixed = NDSolve[Join[Equations, IniConds], {x1, x2}, {t, 0, tMax},
> Method -> CRK4, MaxSteps -> 1000000]
> AbsoluteTime[] - tt0

You could use StartingStepSize -> 0.04 to control the step size.

>
> x1t[t_] := x1[t] /. fixed[];
> x2t[t_] := x2[t] /. fixed[];
>
> ParametricPlot[{x1t[t], x2t[t]}, {t, 0, tMax}]
>
> Out=3.0000000
>
> ***************************************************************************=
> *
>

Just to mention this:

AbsoluteTiming[
var = NDSolve[Join[Equations, IniConds], {x1, x2}, {t, 0, tMax},
Method -> "ExplicitRungeKutta"]
]

is quite efficient and uses adaptive step sizes.

Oliver

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

```

• Prev by Date: Re: Rational[a,b] vs Rational[1,2]
• Next by Date: Re: Vector Runge-Kutta ODE solver with compilation?
• Previous by thread: Re: Vector Runge-Kutta ODE solver with compilation?
• Next by thread: Re: Vector Runge-Kutta ODE solver with compilation?