Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116651] Re: Vector Runge-Kutta ODE solver with compilation?
- From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
- Date: Tue, 22 Feb 2011 06:25:45 -0500 (EST)
On Mon, 21 Feb 2011, DmitryG 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 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. Hello Dmitry, > > 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? The compiler returns tensors. So you have two choices. One is to return a list of this structure {{t0,x0,y0},{t1,x1,y1},....} which should be the first choice. Sometimes, however, you have to revert to another mechanism: Say that t is not real but integer then instead of cf = Compile[{{x, _Real, 0}}, {{1, 2. + x, 3.}, {2, 3. + x, 4.}}] cf[1.] you could do something like ll cf2 = Compile[{{x, _Real, 0}}, ll = {1, 2}; {{2. + x, 3.}, {3. + x, 4.}}]; cf2[1.] {{3., 3.}, {4., 4.}} ll {1,2} Hope this helps, Oliver > > Best, > > Dmitry > >