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