Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116667] Re: Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag at gmail.com>
- Date: Wed, 23 Feb 2011 05:24:57 -0500 (EST)
- References: <201102210919.EAA20606@smc.vnet.net> <ijth3r$mck$1@smc.vnet.net>
On 22 Feb., 04:43, DmitryG <einsch... at gmail.com> wrote: > On 21 Feb., 19:29, DmitryG <einsch... at gmail.com> wrote: > > > > > > > > > > > On 21 Feb., 06:08, Oliver Ruebenkoenig <ruebe... at wolfram.com> wrote: > > > > Dear DmityG, > > > > On Mon, 21 Feb 2011, DmitryG wrote: > > > > Dear All, > > > > > Inspired by the new capability of Mathematica 8 to compile pieces of > > > > user code with its own or external C compiler, I am trying to > > > > implement the simple 4th-order Runge-Kutta solver with fixed step in a > > > > vector form, that is, suitable for solving an arbitrary number of > > > > equations. > > > > > The solver without compilation, the idea of which is taken from > > > >http://www.jasoncantarella.com, is working fine: > > > > > (* Definition *) > > > > RK4[f_,x0_,t0_,tMax_,n_] := Module[{h,K1,K2, K3, > > > > K4,x=x0,SolList={{t0,x0}}}, > > > > > h = (tMax - t0)/n; > > > > > 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=Append[SolList,{t,x}] > > > > > ,{k,1,n}]; > > > > SolList > > > > ]; > > > > > (* Testing for a system of two equations *) > > > > F[t_,x_] := {x[[2]],1 - x[[1]] + 4/(x[[1]]^3)}; > > > > t0=AbsoluteTime[]; > > > > Solution=RK4[ F,{1.0,1.0},0.0,200.0,5000]; > > > > AbsoluteTime[]-t0 > > > > ListLinePlot[Take[Solution[[All,2]],100],PlotMarkers- > > > >> Automatic,PlotStyle->{Red}] > > > > > Out[55]= 0.5468750 > > > > This has the problem that it is not working with PackedArrays which wil l > > > not deliver the results with good performance. > > > > Developer`PackedArrayQ@Solution > > > False > > > > Also, note that you store the time but remove it with Solution[[All,2]]. > > > > > The output above is execution time. The original code uses NestList > > > > instead of my Do cycle but the execution time is the same. > > > > > Unfortunately, the compiled version of this code does not work as > > > > expected: > > > > It works as expected, but, for example, the input f needs to be a real > > > vector - not a function definition. There is a difference if you have a > > > function that returns a vector or a (compiled) function that expects a > > > vector as input. > > > > > (* Definition *) > > > > RK4Comp = > > > > Compile[{{f, _Real, 1}, {x0, _Real, 1}, {t0, _Real}, {tMax, _Real}, > > > > {n, _Integer}}, > > > > Module[{h, K1, K2, K3, K4, SolList = {{t0, x0}}, x = x0, t}, > > > > > h = (tMax - t0)/n; > > > > > 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 = Append[SolList, {t, x}] > > > > , {k, 1, n}]; > > > > SolList > > > > ](*,CompilationTarget->"C"*) > > > > ] > > > > > (* Testing for a system of two ODEs *) > > > > F[t_, x_] := {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)}; > > > > t0 = AbsoluteTime[]; > > > > Solution = RK4Comp[ F, {1.0, 1.0}, 0.0, 200.0, 5000]; > > > > AbsoluteTime[] - t0 > > > > ListLinePlot[Take[Solution[[All, 2]], 100], PlotMarkers -> Automatic, > > > > PlotStyle -> {Red}] > > > > > During evaluation of In[57]:= CompiledFunction::cfta: Argument F at > > > > position 1 should be a rank 1 tensor of machine-size real numbers. > > > > That's what this message states. The fact the this still runs is because > > > the compiled code will always revert to the original code when the > > > compiled code can not evaluate. > > > > > Out[61]= 0.5312500 > > > > > Mathematica complains and seemingly proceeds without compilation > > > > because execution time is the same. > > > > > Anybody has an idea of what is happening and how it can be corected? I > > > > believe the problem should be relevant for many. > > > > Here is how I'd do it: > > > > Since this is Mathematica, I'd write a function that will generate me the > > > compiled function like so: > > > > 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[{0., 0.}, {n + 1}]; > > > SolList[[1]] = x0; > > > 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]] = x; > > > , {k, 2, n + 1}]; > > > SolList], CompilationTarget -> "C"]]; > > > > Note, that I have thrown away the input f and instead inject the code= f= > or > > > the ode into the compiled function. (That's the With[] stuff). Also, = I = > do > > > not store t - if you really need that (perhaps because you'd like to > > > implement an adaptive RK) I can show you how to do that too. Also, I > > > replaced the AppendTo with a Table and set the parts of it directly. > > > (Small note: ConstantArray will not to the trick as it does not compi= le= > a= > > s > > > you could see with CompilePrint) > > > > OK, then > > > > f = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)}]; > > > RK4Comp = makeCompRK[f]; > > > > t0 = AbsoluteTime[]; > > > Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 500000]; > > > AbsoluteTime[] - t0 > > > > runs in about 0.75 sec. on my laptop - that is for n = 500000 (and = no= > t > > > n=5000) > > > > Developer`PackedArrayQ@Sol > > > True > > > > Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 5000]; > > > ListLinePlot[Take[Sol, 100], PlotMarkers -> Automatic, > > > PlotStyle -> {Red}] > > > > I hope this helps, > > > > Oliver > > > > > Thank you for your time, > > > > > Dmitry > > > Dear Oliver, > > > thank you for your help! Your code works on my computer. > > > I see there is a whole bunch of problems and I had absolutely no > > chance to solve any of them myself! Mathematica documentation in Help > > and online gives only simple examples of compilations where > > compilation is probably not necessary. Vector RK4 is an important > > generic real-life problem and, I think, its implementation has to be > > included in Mathematica documentation. > > > I see the idea to compile RK4 with the particular system of ODEs. I > > was trying to compile a RK4 that takes any system of ODEs but it did > > not work. I believe, in any case, RK4 compiled with a particular > > system of ODEs should be faster. > > > Of course, I do need t in the solution. My attempt to return to my > > method of building the solution list with SolList = Append[SolList, > > {t, x}] does not work. It seems compiler does not understand lists of > > a more complicated structure. Also modification of your method of > > initiatingSolList to zeros and then defining SolList[[k]]={t, x} does > > not work. In fact, whatever I do to put t into SolList, it does not > > compile. Could you, help with this, please? > > > Best, > > > Dmitry > > Well, after many tries I have written a compilable vector RK4 ODE > solver that stored time t in the output list: > > *********************************************************************** > (* 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[Prepend[x0, t0], {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]] = Prepend[x, t ]; (* This line did not w= ork for > some reason, be careful! *) > , {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[[All, {2, 3}]], nTake], > PlotMarkers -> Automatic, PlotStyle -> {Red}] > > ListLinePlot[{Take[Sol[[All, {1, 2}]], nTake], > Take[Sol[[All, {1, 3}]], nTake]}, PlotMarkers -> Automatic, > PlotStyle -> {Blue, Green}] > > *************************************************************************= ** > > My idea was to extend the output list SolList as {x[[1]],x[[2]],...} - > > > {t,x[[1]],x[[2],...}. I had a lot of problems with the command > > SolList[[k + 1]] = Prepend[x, t ] > > that gave errors and non-compilation for no apparent reason but then > suddenly worked. This format of data storage is convenient to extract > individual curves x[[i]][t] as plotting shows. > > However, storing t leads to doubling of the execution time in > comparison to the method without saving t. This means, probably, that > > SolList[[k + 1]] = Prepend[x, t ] > > consumes so much time as everything else! So I would recommend not to > store the time that can be generated after the problem has been > solved. > > Best, > > Dmitry > > P.S.: A problem remains for me why the same method used as a plugin to > NDSolve is so much slower and less acurate. 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]; (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.
- Follow-Ups:
- Re: Re: Vector Runge-Kutta ODE solver with compilation?
- From: Oliver Ruebenkoenig <ruebenko@wolfram.com>
- Re: Re: Vector Runge-Kutta ODE solver with compilation?
- References:
- Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag@gmail.com>
- Vector Runge-Kutta ODE solver with compilation?