Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116786] Re: Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag at gmail.com>
- Date: Mon, 28 Feb 2011 04:59:43 -0500 (EST)
- References: <ikd6c8$pvr$1@smc.vnet.net>
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