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

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

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