Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116650] Re: Vector Runge-Kutta ODE solver with compilation?
- From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
- Date: Tue, 22 Feb 2011 06:25:34 -0500 (EST)
On Mon, 21 Feb 2011, DmitryG wrote: > 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 > 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 You could use StartingStepSize -> 0.04 to control the step size. > > x1t[t_] := x1[t] /. fixed[[1]]; > x2t[t_] := x2[t] /. fixed[[1]]; > > ParametricPlot[{x1t[t], x2t[t]}, {t, 0, tMax}] > > Out[61]=3.0000000 > > ***************************************************************************= > * > Just to mention this: AbsoluteTiming[ var = NDSolve[Join[Equations, IniConds], {x1, x2}, {t, 0, tMax}, Method -> "ExplicitRungeKutta"] ] is quite efficient and uses adaptive step sizes. Oliver > 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 > >