MathGroup Archive 2011

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

Search the Archive

Re: Vector Runge-Kutta ODE solver with compilation?

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

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[[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 a=
s
> 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

Dear Oliver,

thank you for your help! Your code works on my computer.

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?

Best,

Dmitry


  • Prev by Date: Re: Vector Runge-Kutta ODE solver with compilation?
  • Next by Date: Re: Delete elements from list..
  • Previous by thread: Re: Vector Runge-Kutta ODE solver with compilation?
  • Next by thread: Re: Vector Runge-Kutta ODE solver with compilation?