MathGroup Archive 2011

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

Search the Archive

Re: Vector Runge-Kutta ODE solver with compilation?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg116887] Re: Vector Runge-Kutta ODE solver with compilation?
  • From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
  • Date: Thu, 3 Mar 2011 05:59:14 -0500 (EST)

Dmitry,

On Wed, 2 Mar 2011, DmitryG wrote:

> Continuing efforts to vectorize NDSolve. Now the vectorized code
>
> *********************************************************************
> NN = 1000;  tMax = 50;
> x0 = Table[RandomReal[{0, 1}], {i, 1, NN}];
> Equations = x'[t] == - x[t]/(1 + 300 Total[x[t]]^2/NN^2);

Try to use NN = 10

and then look at Equations.

A way to make that work is for example:

NN = 1000; tMax = 50;
f[x_List] := 1 + 300*Total[x]^2/NN^2
Equations = x'[t] == -x[t]/f[x[t]]

x0 = Table[RandomReal[{0, 1}], {i, 1, NN}];
Timing[Solution =
   NDSolve[{Equations, x[0] == x0}, x, {t, 0, tMax},
    MaxSteps -> 1000000, SolveDelayed -> True]]

xt[t_] := Evaluate[x[t] /. Solution[[1]]];

Plot[{xt[t][[1]], xt[t][[2]], xt[t][[3]]}, {t, 0, 50},
  PlotStyle -> {{Thick, Red}, {Thick, Green}, {Thick, Blue}},
  PlotRange -> {0, 1}]

This seems to be slower than the NDSolve code you give below.

Hope this helps,

Oliver

>
> Timing[Solution = NDSolve[{Equations, x[0] == x0}, x, {t, 0, tMax},
> MaxSteps -> 1000000]]
>
>
> xt[t_] := x[t] /. Solution[[1]];
>
> Plot[{xt[t][[1]], xt[t][[2]], xt[t][[3]]}, {t, 0, 50}, PlotStyle ->
> {{Thick, Red}, {Thick, Green}, {Thick, Blue}}, PlotRange -> {0, 1}]
>
> *****************************************************************************
> computes something, the solution x[t] is a vector, but the resulting
> plot differs from the plot generated by the non-vectorized version
>
> ******************************************************************************
> NN = 1000;  tMax = 50;
> x0 = Table[RandomReal[{0, 1}], {i, 1, NN}];
> IniConds = Table[x[i][0] == x0[[i]], {i, 1, NN}];
> Vars = Table[x[i], {i, 1, NN}];
> Timing[Equations = Table[x[i]'[t] == -x[i][t]/(1 + 300 Sum[x[j][t],
> {j, 1, NN}]^2/NN^2), {i, 1,NN}];]
>
> Timing[Solution = NDSolve[Join[Equations, IniConds], Vars, {t, 0,
> tMax}, MaxSteps -> 1000000)];
>
> x1t[t_] := x[1][t] /. Solution[[1]];
> x2t[t_] := x[2][t] /. Solution[[1]];
> x3t[t_] := x[3][t] /. Solution[[1]];
> Plot[{x1t[t], x2t[t], x3t[t]}, {t, 0, tMax}, PlotStyle -> {{Thick,
> Red}, {Thick, Green}, {Thick, Blue}}]
>
> **************************************************************************
>
> and the RK4 routine above. Although the vectorized program is much
> faster, its results are wrong. Where is the error??
>
> Best,
>
> Dmitry
>
>


  • Prev by Date: TransformedDistribution
  • Next by Date: Re: Problem exporting to LaTex
  • Previous by thread: Re: Vector Runge-Kutta ODE solver with compilation?
  • Next by thread: Re: Vector Runge-Kutta ODE solver with compilation?