Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116706] Re: Vector Runge-Kutta ODE solver with compilation?
- From: Leonid Shifrin <lshifr at gmail.com>
- Date: Thu, 24 Feb 2011 06:24:41 -0500 (EST)
Hi Oliver,
Just a side remark: if using delayed assignment is acceptable for RHS,
another way to avoid
error messages (and generally the evaluation of RHS) would be
In[48]:=
RHS:={x[[2]],1-x[[1]]+4/x[[1]]^3};
Function@@Prepend[Hold[RHS]/.OwnValues[RHS],Unevaluated[{t,x}]]
Out[49]= Function[{t,x},{x[[2]],1-x[[1]]+4/x[[1]]^3}]
I found this idiom useful on several occasions.
Cheers,
Leonid
On Wed, Feb 23, 2011 at 2:25 PM, Oliver Ruebenkoenig
<ruebenko at wolfram.com>wrote:
>
> On Wed, 23 Feb 2011, DmitryG wrote:
>
> [....]
>
>
> > Dear Oliver and All interested,
> >
> > Today I tried to apply the above very efficient vector RK4 method with
> > compilation to my real problems that I was solving with NDSolve
> > before. Unfortunately, it has proven to be more complicated than it
> > seemed. To illustrate the difficulties, I will use the same working
> > model as above and will try to generalize it. So, the working code is
> >
> > *************************************************
> >
> > (* Definition of the RK4 solver *)
> > makeCompRK[fIn_] :=
> > With[{f = fIn},
> > Compile[{{x0, _Real, 1}, {t0, _Real}, {tMax, _Real}, {n,
> > _Integer}},
> > Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
> > h = (tMax - t0)/n;
> > SolList = Table[x0, {n + 1}];
> > 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[[k + 1]] = x
> >
> > , {k, 1, n}];
> > SolList], CompilationTarget -> "C"]];
> >
> > (* Calculation *)
> > F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
> > RK4Comp = makeCompRK[F];
> > x0 = {1.0, 1.0}; t0 = 0; tMax = 200; n = 5000;
> > tt0 = AbsoluteTime[];
> > Sol = RK4Comp[x0, t0, tMax, n];
> > AbsoluteTime[] - tt0
> > Developer`PackedArrayQ@Sol
> >
> > (* Plotting *)
> > nTake = 200;
> > ListLinePlot[Take[Sol, nTake], PlotMarkers -> Automatic, PlotStyle ->
> > {Red}] (* Trajectory *)
> >
> > (* Time dependences *)
> > x1List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 1]]}, {k, 0, n}];
> > x2List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 2]]}, {k, 0, n}];
> > ListLinePlot[{Take[x1List, nTake], Take[x2List, nTake]},
> > PlotMarkers -> Automatic, PlotStyle -> {Blue, Green}]
> >
> ***************************************************************************=
> > **************************
> >
> > In this example, the vector of ODE's RHS' is defined by
> >
> > F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
> >
> > In real life, this is something more complicated, a result of a many-
> > step calculation. One can try to generalize the code as
> >
> > RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3};
> > F = Function[{t, x}, RHS];
>
> Dmitry,
>
> your friend is - in princial -
>
> F = Function[{t, x}, Evaluate[RHS]]
>
> That has a small problem that we used Part in the RHS and Part gives a
> message that you try to access elements that the _symbol_ x does not have.
>
> One way to work around this is
>
> RHS = Quiet[{x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
> F = Function[{t, x}, Evaluate[RHS]];
>
>
> If you evaluate F you'll see that the RHS in now in the function body. If
> you look at the F you posted, there the body is RHS, which is not what
> you want. All this hast to do with
>
> Attributes[Function]
>
> HoldAll.
>
> Oliver
>
>
>
> >
> > (that should be equivalent to the above) and then substitute RHS by
> > whatever you have for your problem. However, Mathematica complains and
> > produces crap results. My suspicion is that it is only F =
> > Function[{t, x}, RHS]; that is injected into the compiled RK4 code
> > whereas RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; is not injected.
> > However, I cannot be sure because it was never explained.
> >
> > Further idea is to define F as a procedure and put everything inside
> > it. A toy example is
> >
> > F = Function[{t, x}, RHS = {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; RHS];
> >
> > This does not compile, although then retreating to the non-compilation
> > mode produces correct results. My experiments show that you can only
> > use a single expression in the body of F that does not use any
> > definitions. For instance,
> >
> > F = Function[{t, x}, Table[-x[[i]]^i, {i, 1, 2}]];
> >
> > compiles and works fine but this is, again, no more than a toy
> > example.
> >
> > To summarize, I am very excited by the prospect to speed up my
> > computations by a factor about 100 but I don't see how to generalize
> > the approach for real-life problems.
> >
> >
> >
>
>