Re: Vector Runge-Kutta ODE solver with compilation?

*To*: mathgroup at smc.vnet.net*Subject*: [mg116706] Re: Vector Runge-Kutta ODE solver with compilation?*From*: Leonid Shifrin <lshifr at gmail.com>*Date*: Thu, 24 Feb 2011 06:24:41 -0500 (EST)

Hi Oliver, Just a side remark: if using delayed assignment is acceptable for RHS, another way to avoid error messages (and generally the evaluation of RHS) would be In[48]:= RHS:={x[[2]],1-x[[1]]+4/x[[1]]^3}; Function@@Prepend[Hold[RHS]/.OwnValues[RHS],Unevaluated[{t,x}]] Out[49]= Function[{t,x},{x[[2]],1-x[[1]]+4/x[[1]]^3}] I found this idiom useful on several occasions. Cheers, Leonid On Wed, Feb 23, 2011 at 2:25 PM, Oliver Ruebenkoenig <ruebenko 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. > > > > > > > >