Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116869] Re: Vector Runge-Kutta ODE solver with compilation?
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Thu, 3 Mar 2011 05:55:58 -0500 (EST)
Mathematica doesn't understand that x[t] is a vector, because it isn't --
not in that code, at least. It's not even defined.
To define it from the NDSolve solution:
NN=100;tMax=50;
x0=Table[RandomReal[{0,1}],{i,1,NN}];
ClearAll[x]
equations=x'[t]==-x[t]/(1+300 Total[x[t]]^2/NN^2);
Timing[solution=NDSolve[{equations,x[0]==x0},x,{t,0,tMax},MaxSteps->1000000]]
x[t_?NumericQ]=x[t]/.First@solution;
{0.006209,{{x->InterpolatingFunction[{{0.,50.}},<>]}}}
x1t[t_] := x[t][[1]]
x2t[t_] := x[t][[2]]
x3t[t_] := x[t][[3]]
Plot[{x1t[t], x2t[t], x3t[t]}, {t, 0, tMax},
PlotStyle -> {Red, Green, Blue}]
or
Clear[p]
colors = {Red, Green, Blue};
p[i_] := Plot[x[t][[i]], {t, 0, tMax}, PlotStyle -> colors[[i]]]
Show[Array[p, 3]]
or
Clear[p]
colors = {Red, Green, Blue};
p[i_] := Plot[x[t][[i]], {t, 0, tMax},
PlotStyle -> colors[[Mod[i, 3, 1]]]]
Show[Array[p, Length@x0]]
or better yet,
Clear[p]
color[i_] := ColorData["BrightBands"][x0[[i]]]
p[i_] := Plot[x[t][[i]], {t, 0, tMax}, PlotStyle -> color@i]
Show[Array[p, Length@x0]]
You COULD make most of the original code work, this way:
NN=100;tMax=50;
x0=Table[RandomReal[{0,1}],{i,1,NN}];
Clear[x]
equations=x'[t]==-x[t]/(1+300 Total[x[t]]^2/NN^2);
Timing[solution=NDSolve[{equations,x[0]==x0},x,{t,0,tMax},MaxSteps->1000000]]
{0.006306,{{x->InterpolatingFunction[{{0.,50.}},<>]}}}
x1t[t_]:=(x[t]/.solution[[1]])[[1]]
x2t[t_]:=(x[t]/.solution[[1]])[[2]]
x3t[t_]:=(x[t]/.solution[[1]])[[3]]
Plot[{x1t[t],x2t[t],x3t[t]},{t,0,tMax},PlotStyle->{Red,Green,Blue}]
Bobby
On Wed, 02 Mar 2011 03:33:02 -0600, DmitryG <einschlag at gmail.com> wrote:
> Dear Oliver and Daniel,
>
> Vectorizing RHS is a magic trick that did the job, thank you a lot!
> Now the compiled code is short and its length does not depend on NN.
> The execution time strongly decreased, as one can expect from vector
> operations. I can also compile in C without problems. Strangely, the
> execution time is longer in the case of compilation in C!
>
> I hope I will be able to do vectorization for my real problem where
> RHS is similar but more complicated.
>
> Another idea is to do vectorization for NDSolve. I have always
> generated equations as Lists such as Equations=Table[...] but maybe I
> could completely vectorize NDSolve and achieve a great speed? I have
> tried with the system of equations we are discussing
>
> ******************************************************************************
>
> NN = 100; 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[0] == x0}, x, {t, 0, tMax},
> MaxSteps -> 1000000]]
>
> x1t[t_] := x[t][[1]] /. Solution[[1]];
> x2t[t_] := x[t][[2]] /. Solution[[1]];
> x3t[t_] := x[t][[3]] /. Solution[[1]];
> Plot[{x1t[t], x2t[t], x3t[t]}, {t, 0, tMax}, PlotStyle -> {Red, Green,
> Blue}]
>
> **********************************************************************************
>
> but this is not working. Probably Mathematica does not understand that
> x[t] is a vector? How to fix the problem?
>
> Thank you for your help,
>
> Dmitry
>
>
--
DrMajorBob at yahoo.com