Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Vector Runge-Kutta ODE solver with compilation?

• To: mathgroup at smc.vnet.net
• Subject: [mg116631] Re: Vector Runge-Kutta ODE solver with compilation?
• From: DmitryG <einschlag at gmail.com>
• Date: Tue, 22 Feb 2011 04:43:15 -0500 (EST)
• References: <201102210919.EAA20606@smc.vnet.net> <ijth3r$mck$1@smc.vnet.net> <ijv00n$660$1@smc.vnet.net>

On 21 Feb., 19:29, DmitryG <einsch... at gmail.com> 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 wil=
l
> > 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]]=

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