[Date Index]
[Thread Index]
[Author Index]
Re: Vector Runge-Kutta ODE solver with compilation?
*To*: mathgroup at smc.vnet.net
*Subject*: [mg116611] Re: Vector Runge-Kutta ODE solver with compilation?
*From*: DmitryG <einschlag at gmail.com>
*Date*: Mon, 21 Feb 2011 19:29:12 -0500 (EST)
*References*: <201102210919.EAA20606@smc.vnet.net> <ijth45$mcp$1@smc.vnet.net>
On 21 Feb., 06:08, Oliver Ruebenkoenig <ruebe... at wolfram.com> wrote:
> 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
Oliver,
This looks nice and works but is very slow. What is the step size h in
the fixed-step RK4 algorithm? Here is my working version for the same
equations:
**************************************************************************
CRK4[]["Step"[rhs_, t_, h_, y_, yp_]] := Module[{k0, k1, k2, k3},
k0 = h yp;
k1 = h rhs[t + h/2, y + k0/2];
k2 = h rhs[t + h/2, y + k1/2];
k3 = h rhs[t + h, y + k2];
{h, (k0 + 2 k1 + 2 k2 + k3)/6}]
CRK4[___]["DifferenceOrder"] := 4
CRK4[___]["StepMode"] := Fixed (* What is "Fixed"?? Mathematica does
not know this. What is the step h?? *)
Equations = {x1'[t] == x2[t], x2'[t] == 1 - x1[t] + 4/(x1[t]^3)};
IniConds = {x1[0] == 1, x2[0] == 1};
tMax = 200;
tt0 = AbsoluteTime[];
fixed = NDSolve[Join[Equations, IniConds], {x1, x2}, {t, 0, tMax},
Method -> CRK4, MaxSteps -> 1000000]
AbsoluteTime[] - tt0
x1t[t_] := x1[t] /. fixed[[1]];
x2t[t_] := x2[t] /. fixed[[1]];
ParametricPlot[{x1t[t], x2t[t]}, {t, 0, tMax}]
Out[61]=3.0000000
***************************************************************************=
*
Here, for tMax=500, the execution time is 3.00 and the plot shows that
the precision is insufficient because of a too large h. For larger
tMax it becomes dramatic.
Best,
Dmitry
Prev by Date:
**Re: Chain of packets in MathLink: are packets always strictly ordered?**
Next by Date:
**Re: Vector Runge-Kutta ODE solver with compilation?**
Previous by thread:
**Re: Vector Runge-Kutta ODE solver with compilation?**
Next by thread:
**Re: Vector Runge-Kutta ODE solver with compilation?**
| |