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: [mg116602] Re: Vector Runge-Kutta ODE solver with compilation?
• From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
• Date: Mon, 21 Feb 2011 06:08:08 -0500 (EST)
• References: <201102210919.EAA20606@smc.vnet.net>

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

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

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


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[[1]] = 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 as 
you could see with CompilePrint)

OK, then

f = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/(x[[1]]^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
>
>

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