Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116898] Re: Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag at gmail.com>
- Date: Fri, 4 Mar 2011 03:37:01 -0500 (EST)
- References: <iknsbs$k8o$1@smc.vnet.net>
On 3 Mrz., 06:00, Oliver Ruebenkoenig <ruebe... at wolfram.com> wrote:
> 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
Thank you a lot, Oliver!
Why have you added the option SolveDelayed->True in NDSolve? Is it
important here or in general? In the help it is written this option is
useful in the case of singularities, to let Mathematica do some
algebra.
Best,
Dmitry