MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Vector Runge-Kutta ODE solver with compilation?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg116741] Re: Vector Runge-Kutta ODE solver with compilation?
  • From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
  • Date: Fri, 25 Feb 2011 06:38:09 -0500 (EST)

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 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.
>
> 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."

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
>
>


  • Prev by Date: Re: FinancialData Function Not Working for Property "Members"
  • 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?