Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116822] Re: Vector Runge-Kutta ODE solver with compilation?
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Tue, 1 Mar 2011 05:23:48 -0500 (EST)
DmitryG wrote:
> On 27 Feb., 04:43, Oliver Ruebenkoenig <ruebe... at wolfram.com> wrote:
>
>> You could use CompilePrint to see that for NN=11 this generates a some 900
>> line so code - quite a lot of which are Part[]. Make sure you understand
>> what you are doing here - you in-line this sum into the code.
>>
>>
> ......................
>> Weird code? - No, not really.
>> ...
>>
>
> Hi Oliver,
>
> In this short program I mimick my real-life program that solves a
> system of NN ODE's with the maximal coupling: the RHS of each of NN
> equations contains all NN variables. Any molecular dynamics (MD)
> simulation is of this kind: Acceleration of any particle depends on
> positions of all particles.
>
> This application of Mathematica's compilation is crucial but there are
> two problems that I have encountered:
>
> 1) Compilation in C is prohibitively slow above NN=50 or so. In this
> case Mathematica's own compiler is still doing a good job.
>
> 2) For NN=1000 or above Mathematica's compiler consumes all RAM (my
> laptop has 3 GB RAM) and quits. This is very disappointing because
> generating the list of RHS does not take much memory and the RK 4
> routine does not consume much memory as well. Thus it would be only a
> matter of computer time to solve a problem for any NN. However,
> compiler generates a blown-up code (component by component, with the
> vector structure abandoned) that does not fit into memory.
>
> Thank you for suggestion to look into the compiled code with
> CompilePrint. Yes, the code is very long, and its structure is the
> same for the initial problem we considered, the problem with
>
> RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]
>
> Also in this case the compiled code contains Part[].
>
> Is there any way to optimize this calculation or we have an intrinsic
> limitation of Mathematica here? Below is the code again:
>
> *******************************************************************************
> (* 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=15;
> Timing[RHS=Quiet[Table[-x[[i]]/(1+300Sum[x[[j]],{j,1,NN}]^2/NN^2),{i,
> 1,NN}]];] (* Quiet to suppress complaints *)
> F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
> RHS *)
> Print[]
>
> tt0=AbsoluteTime[];
> Timing[RK4Comp=makeCompRK[F];]
> AbsoluteTime[]-tt0
> Print[];
>
> x0=Table[RandomReal[{0,1}],{i,1,NN}]; t0=0; tMax=50; 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}]
>
>
> *********************************
>
> Thank you for your help,
>
> Dmitry
>
I'll show code below that seems to be better behaved. I believe it is
functionally the same as yours but you'll need to check that. Itt is
based on some things Oliver Ruebenkoenig suggested (possibly you had not
yet seen that post at the time you wrote this). The main change is to
vectorize that RHS. I also use a Table instead of loop to fill in the
solution, but that probably is not too important.
makeCompRK[fIn_] :=
With[{f = fIn},
Compile[{{x0, _Real,
1}, {t0, _Real}, {tMax, _Real}, {n, _Integer}},
Module[{h, K1, K2, K3, K4, x = x0, t},
h = (tMax - t0)/n;
Prepend[Table[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,
{k, 1, n}], x0]
](*,CompilationTarget->"C"*)]];
(*Calculation*)
NN = 150;
Timing[RHS =
Quiet[-x/(1 +
300 Total[x]^2/NN^2)];] (*Quiet to suppress complaints*)
F =
Function[{t, x}, Evaluate[RHS]];(*Evaluate to use actual form of RHS*)
Out[136]= {0., Null}
tt0 = AbsoluteTime[];
Timing[RK4Comp = makeCompRK[F];]
AbsoluteTime[] - tt0
x0 =
Table[RandomReal[{0, 1}], {i, 1, NN}]; t0 = 0; tMax = 50; n = 500;
tt0 = AbsoluteTime[];
Sol = RK4Comp[x0, t0, tMax, n];
AbsoluteTime[] - tt0
Out[144]= 0.004684
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}]
With this version I had no difficulty handling NN=1500 either using C or
the virtual machine interpreter of Compile.
Daniel Lichtblau
Wolfram Research