       Re: Vector Runge-Kutta ODE solver with compilation?

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

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

> On 21 Feb., 06:08, Oliver Ruebenkoenig <ruebe... at wolfram.com> wrote:
>> Dear DmityG,
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> 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
>>
>> This has the problem that it is not working with PackedArrays which will
>> not deliver the results with good performance.
>>
>> Developer`PackedArrayQ@Solution
>> False
>>
>> Also, note that you store the time but remove it with Solution[[All,2]].
>>
>>
>>
>>> 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:
>>
>> It works as expected, but, for example, the input f needs to be a real
>> vector - not a function definition. There is a difference if you have a
>> function that returns a vector or a (compiled) function that expects a
>> vector as input.
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>> (* 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. >>
>>
>> That's what this message states. The fact the this still runs is because
>> the compiled code will always revert to the original code when the
>> compiled code can not evaluate.
>>
>>> 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.
>>
>> Here is how I'd do it:
>>
>> Since this is Mathematica, I'd write a function that will generate me the
>> compiled function like so:
>>
>> makeCompRK[fIn_] := With[{f = fIn},
>>     Compile[{{x0, _Real,
>>        1}, {t0, _Real}, {tMax, _Real}, {n, _Integer}},
>>      Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
>>       h = (tMax - t0)/n;
>>       SolList = Table[{0., 0.}, {n + 1}];
>>       SolList[] = x0;
>>       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[[k]] = x;
>>        , {k, 2, n + 1}];
>>       SolList], CompilationTarget -> "C"]];
>>
>> Note, that I have thrown away the input f and instead inject the code for
>> the ode into the compiled function. (That's the With[] stuff). Also, I do
>> not store t - if you really need that (perhaps because you'd like to
>> implement an adaptive RK) I can show you how to do that too. Also, I
>> replaced the AppendTo with a Table and set the parts of it directly.
>> (Small note: ConstantArray will not to the trick as it does not compile a=
> s
>> you could see with CompilePrint)
>>
>> OK, then
>>
>> f = Function[{t, x}, {x[], 1 - x[] + 4/(x[]^3)}];
>> RK4Comp = makeCompRK[f];
>>
>> t0 = AbsoluteTime[];
>> Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 500000];
>> AbsoluteTime[] - t0
>>
>> runs in about 0.75 sec. on my laptop - that is for n = 500000 (and not
>> n=5000)
>>
>> Developer`PackedArrayQ@Sol
>> True
>>
>> Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 5000];
>> ListLinePlot[Take[Sol, 100], PlotMarkers -> Automatic,
>>   PlotStyle -> {Red}]
>>
>> I hope this helps,
>>
>> Oliver
>>
>>
>>
>>
>>
>>
>>
>>> Thank you for your time,
>>
>>> Dmitry
>
> Dear Oliver,
>
> thank you for your help! Your code works on my computer.

Hello Dmitry,

>
> I see there is a whole bunch of problems and I had absolutely no
> chance to solve any of them myself! Mathematica documentation in Help
> and online gives only simple examples of compilations where
> compilation is probably not necessary. Vector RK4 is an important
> generic real-life problem and, I think, its implementation has to be
> included in Mathematica documentation.
>
> I see the idea to compile RK4 with the particular system of ODEs. I
> was trying to compile a RK4 that takes any system of ODEs but it did
> not work. I believe, in any case, RK4 compiled with a particular
> system of ODEs should be faster.
>
> Of course, I do need t in the solution. My attempt to return to my
> method of building the solution list with SolList = Append[SolList,
> {t, x}] does not work. It seems compiler does not understand lists of
> a more complicated structure. Also modification of your method of
> initiatingSolList to zeros and then defining SolList[[k]]={t, x} does
> not work. In fact, whatever I do to put t into SolList, it does not
> compile. Could you, help with this, please?

The compiler returns tensors. So you have two choices. One is to return a
list of this structure

{{t0,x0,y0},{t1,x1,y1},....} which should be the first choice.

Sometimes, however, you have to revert to another mechanism: Say that t is
not real but integer then instead of

cf = Compile[{{x, _Real, 0}}, {{1, 2. + x, 3.}, {2, 3. + x, 4.}}]
cf[1.]

you could do something like

ll

cf2 = Compile[{{x, _Real, 0}},
ll = {1, 2}; {{2. + x, 3.}, {3. + x, 4.}}];

cf2[1.]
{{3., 3.}, {4., 4.}}

ll
{1,2}

Hope this helps,
Oliver

>
> Best,
>
> Dmitry
>
>

```

• Prev by Date: Re: Vector Runge-Kutta ODE solver with compilation?
• Next by Date: Re: weibull plot on weibull scaled paper
• Previous by thread: Re: Vector Runge-Kutta ODE solver with compilation?
• Next by thread: Re: Vector Runge-Kutta ODE solver with compilation?