Services & ResourcesWolfram ForumsMathGroup Archive

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

Vector Runge-Kutta ODE solver with compilation?

• To: mathgroup at smc.vnet.net
• Subject: [mg116585] Vector Runge-Kutta ODE solver with compilation?
• From: DmitryG <einschlag at gmail.com>
• Date: Mon, 21 Feb 2011 04:19:25 -0500 (EST)

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

• Follow-Ups:
  • Re: Vector Runge-Kutta ODE solver with compilation?
    • From: Oliver Ruebenkoenig <ruebenko@wolfram.com>
  • Re: Vector Runge-Kutta ODE solver with compilation?
    • From: Oliver Ruebenkoenig <ruebenko@wolfram.com>