Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2011

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

Search the Archive

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


  • Prev by Date: Re: Integral with singularities
  • Next by Date: Re: How to deal with big matrix 2
  • Previous by thread: Re: Vector Runge-Kutta ODE solver with compilation?
  • Next by thread: Re: Vector Runge-Kutta ODE solver with compilation?