Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116715] Re: Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag at gmail.com>
- Date: Fri, 25 Feb 2011 06:33:27 -0500 (EST)
- References: <201102231024.FAA09857@smc.vnet.net> <ik2qs3$avl$1@smc.vnet.net> <ik5f20$sd2$1@smc.vnet.net>
On 24 Feb., 06:22, DmitryG <einsch... at gmail.com> wrote:
> On 23 Feb., 06:25, Oliver Ruebenkoenig <ruebe... 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-compilatio=
n
> > > 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.
>
> Thank you, Oliver, this was a great help.
>
> Now, I am trying to introduce sampling of the results to reduce the
> number of output points. This is also very important. However, there
> is a new problem. Here is the code:
>
> *************************************************************************=
** ***
> (* Vector RK4 routine *)
> makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
> {ttMax,_Real},{n,_Integer},{nsamp,_Integer}},
> Module[{h,K1,K2,K3,K4,SolList,x=x0,t,ksamp},
> h=(ttMax-t0)/n;
> SolList=Table[x0,{nsamp+1}];
> Do[t=t0+k h; ksamp=k nsamp/n;
> 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;
>
> If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>
> ,{k,1,n}];
> SolList],CompilationTarget->"C"]];
>
> (* Calculation *)
> RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]; (* Quiet to suppress
> complaints *)
> F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
> RHS *)
> RK4Comp=makeCompRK[F];
> x0={1.0,1.0}; t0=0; tMax; NSteps=5000; NSampling0=
;
> tt0=AbsoluteTime[];
> Sol=RK4Comp[x0,t0,tMax,NSteps,NSampling];
> AbsoluteTime[]-tt0
> Developer`PackedArrayQ@Sol
> *************************************************************************=
** ***
>
> The problem is in the line
>
> If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>
> In the first execution, Mathematica complains
>
> Compile::cpintlt: "ksamp$+1 at position 2 of SolList$[[ksamp$+1]]
> should be either a nonzero integer or a vector of nonzero integers;
> evaluation will use the uncompiled function."
>
> I do not understand this. Still Mathematica compiles and the results
> are correct.
>
> In the second execution, Mathematica complains more and does not
> compile, although the results are still correct.
>
> If I use the command
>
> If[IntegerQ[kSampling],SolList[[kSampling+1]]=x
>
> Mathematica complains again:
>
> CompiledFunction::cfse: "Compiled expression False should be a machine-
> size real number)."
>
> does not compile, but the results are correct.
>
> What is wrong here? I believe I am doing everything in a natural way.
> But compiling is too tricky..
>
> Dmitry
The command that works is
If[FractionalPart[ksamp] == 0, SolList[[IntegerPart[ksamp] + 1]] x];
(* Although assignment is done for ksamp integer, IntegerPart is
needed! *)
Commands including the test
IntegerQ[ksamp]
do not work.
Dmitry
- References:
- Re: Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag@gmail.com>
- Re: Vector Runge-Kutta ODE solver with compilation?