Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116609] Re: Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag at gmail.com>
- Date: Mon, 21 Feb 2011 19:28:50 -0500 (EST)
- References: <201102210919.EAA20606@smc.vnet.net> <ijth3r$mck$1@smc.vnet.net>
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 will > 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 for > 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 compile 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 not > 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
- References:
- Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag@gmail.com>
- Vector Runge-Kutta ODE solver with compilation?