Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116985] Re: Vector Runge-Kutta ODE solver with compilation?
- From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
- Date: Sun, 6 Mar 2011 05:44:56 -0500 (EST)
On Fri, 4 Mar 2011, DmitryG wrote:
> 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
>
>
>
Sorry, that was a left-over from an experiment. I wanted to try something
a in a little direction but then forgot to remove that option then.
Oliver