Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116749] Re: Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag at gmail.com>
- Date: Sat, 26 Feb 2011 06:06:24 -0500 (EST)
- References: <ik84bo$n4e$1@smc.vnet.net>
On 25 Feb., 06:38, Oliver Ruebenkoenig <ruebe... at wolfram.com> wrote:
> On Thu, 24 Feb 2011, DmitryG 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 h=
ave.
>
> >> 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 an=
d
> >>> 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; NSamplin=
g0;
> > 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."
>
> Of what type is ksamp?
>
> ksamp = k nsamp/n
>
> For example
>
> Head[5 3/2]
>
> In worst case this is a rational and you check that with
> FractionalPart. But Part needs an integer at position 2. With using Floor
> you can tell the compiler that ksamp is always going to be an integer.
>
> SolList[[Floor[ksamp] + 1]]
>
> Oliver
>
>
>
>
>
>
>
>
>
> > 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
Thank you for the response, Oliver! I have already fugured it out.
There is now a problem with compilation in C in the case of a system
of N equations. For N above 10, it takes a long time to compile, and
the compilation time is dramatically increasing with N, so that one
cannot solve any real-life problems with large N using compilation in
C. On the other hand, using Mathematca's own compiler is non-
problematic. In both cases the solution is correct. Here is the code:
***************************************************************************=
***************
(* Definition *)
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 *)
NN=30;
RHS=Quiet[Table[-x[[i]]/(1+Sum[x[[j]],{j,1,NN}]^2),{i,1,NN}]]; (*
Quiet to suppress complaints *)
F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
RHS *)
tt0=AbsoluteTime[];
Timing[RK4Comp=makeCompRK[F];]
AbsoluteTime[]-tt0
Print[];
x0=Table[RandomReal[{0,1}],{i,1,NN}]; t0=0; tMax0; n=500;
tt0=AbsoluteTime[];
Sol=RK4Comp[x0,t0,tMax,n];
AbsoluteTime[]-tt0
Print["Compilation: ",Developer`PackedArrayQ@Sol]
(* Plotting 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}];
x3List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,3]]},{k,0,n}];
ListLinePlot[{x1List,x2List,x3List},PlotMarkers->Automatic,PlotStyle-
>{Blue,Green,Red},PlotRange->{0,1}]
Compile::nogen: A library could not be generated from the compiled
function. >>
Out[4]= {11.172,Null}
Out[5]= 142.0625000
Out[10]= 0.0468750
Compilation: True
***************************************************************************=
**********************************
Whereas execution is fast, 0.0468750, compilation takes a long time. I
measure both Timing and AbsoluteTime. The former, I believe, is the
time used by Mathematica, whereas the latter is the total time used by
Mathematica and my C compiler, Microsoft Visual C++ 2010 Express. Both
times are long. Presumably, Mathematica generates some weird code that
gives a hard time to C compiler. Here I also have an error message
concerning a library but I haven' had this message as I ran this
program before.
Any help?
Dmitry