Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116585] Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag at gmail.com>
- Date: Mon, 21 Feb 2011 04:19:25 -0500 (EST)
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 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: (* 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. >> 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. Thank you for your time, Dmitry
- Follow-Ups:
- Re: Vector Runge-Kutta ODE solver with compilation?
- From: Oliver Ruebenkoenig <ruebenko@wolfram.com>
- Re: Vector Runge-Kutta ODE solver with compilation?
- From: Oliver Ruebenkoenig <ruebenko@wolfram.com>
- Re: Vector Runge-Kutta ODE solver with compilation?