       Re: Vector Runge-Kutta ODE solver with compilation?

• To: mathgroup at smc.vnet.net
• Subject: [mg116882] Re: Vector Runge-Kutta ODE solver with compilation?
• From: DrMajorBob <btreat1 at austin.rr.com>
• Date: Thu, 3 Mar 2011 05:58:19 -0500 (EST)

```The first version is completely wrong, because this is the wrong equation
set:

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)

Derivative[x][t] == -(x[t]/(1 + (3 t^2)/10000))

That's because Total[anyHead[arguments]] is the sum of the arguments. In
this case,

Total[x[t]]

t

Here's a version that works:

NN = 1000; tMax = 50;
x0 = RandomReal[{0, 1}, NN];

Clear[x]
ones = ConstantArray[1, NN];
equations = x'[t] == -x[t]/(1 + 300 (ones.x[t])^2/NN^2);
Timing[solution =
NDSolve[{equations, x == x0}, x, {t, 0, tMax},
MaxSteps -> 1000000];]

{0.048119, Null}

xt[t_] := x[t] /. solution[];
Plot[{xt[t][], xt[t][], xt[t][]}, {t, 0, 50},
PlotStyle -> {{Thick, Red}, {Thick, Green}, {Thick, Blue}},
PlotRange -> {0, 1}]

That seems to match the following plot:

Clear[x]
initial = Table[x[i] == x0[[i]], {i, 1, NN}];
vars = Array[x, 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, initial], vars, {t, 0, tMax},
MaxSteps -> 1000000];
x1t[t_] := x[t] /. solution[];
x2t[t_] := x[t] /. solution[];
x3t[t_] := x[t] /. solution[];
p = Plot[{x1t[t], x2t[t], x3t[t]}, {t, 0, tMax},
PlotStyle -> {{Thick, Red}, {Thick, Green}, {Thick, Blue}}];
]
p

Bobby

On Wed, 02 Mar 2011 03:34:19 -0600, DmitryG <einschlag at gmail.com> 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);
>
> Timing[Solution = NDSolve[{Equations, x == x0}, x, {t, 0, tMax},
> MaxSteps -> 1000000]]
>
>
> xt[t_] := x[t] /. Solution[];
>
> Plot[{xt[t][], xt[t][], xt[t][]}, {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] == 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[t] /. Solution[];
> x2t[t_] := x[t] /. Solution[];
> x3t[t_] := x[t] /. Solution[];
> 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
>

--
DrMajorBob at yahoo.com

```

• Prev by Date: Re: Filtering data from numerical minimization
• Next by Date: Re: Vector Runge-Kutta ODE solver with compilation?
• Previous by thread: Re: Vector Runge-Kutta ODE solver with compilation?
• Next by thread: Re: Vector Runge-Kutta ODE solver with compilation?