Re: Re: Vector Runge-Kutta ODE solver with compilation?

*To*: mathgroup at smc.vnet.net*Subject*: [mg116681] Re: [mg116667] Re: Vector Runge-Kutta ODE solver with compilation?*From*: Oliver Ruebenkoenig <ruebenko at wolfram.com>*Date*: Wed, 23 Feb 2011 06:25:29 -0500 (EST)*References*: <201102231024.FAA09857@smc.vnet.net>

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. > > >

**References**:**Re: Vector Runge-Kutta ODE solver with compilation?***From:*DmitryG <einschlag@gmail.com>