MathGroup Archive 2011

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

Search the Archive

Vector Runge-Kutta ODE solver with compilation?

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

The solver without compilation, the idea of which is taken from, is working fine:

(* Definition *)
RK4[f_,x0_,t0_,tMax_,n_] := Module[{h,K1,K2, K3,

h = (tMax - t0)/n;

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;



(* Testing for a system of two equations *)
F[t_,x_] := {x[[2]],1 - x[[1]] + 4/(x[[1]]^3)};
Solution=RK4[ F,{1.0,1.0},0.0,200.0,5000];

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

(* 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;

     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}];

(* 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,


  • Prev by Date: Re: Color grid with x and y args to visualize effects of 2D
  • Next by Date: Re: weibull plot on weibull scaled paper
  • Previous by thread: Re: Misprint in the Documentation?
  • Next by thread: Re: Vector Runge-Kutta ODE solver with compilation?