Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116928] Re: Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag at gmail.com>
- Date: Fri, 4 Mar 2011 03:42:27 -0500 (EST)
- References: <ikihnp$7sm$1@smc.vnet.net> <ikl30v$sra$1@smc.vnet.net> <iknsjl$kbs$1@smc.vnet.net>
Thank you, Bob, Oliver, and Daniel, for excellent suggestions!
Yes, I have already seen that Total[x[t]] outputs t.
I've tried to tell Mathematica that x[t] is a vector by replacing x[t]
by IdentityMatrix[NN].x[t] but it does not work.
It is interesting that in the version using compiled RK4 routine the
only place that shows x[t] is a vector is the initial condition for
x[t] at t=0, a vector x0. This is sufficient for Mathematica to work
properly. However, in the program using NDSolve setting the vector
initial condition x[0]==x0 is insufficient.
In fact, the real-life equations that I am solving are a little more
complicated and contain AMatr.x[t] in the RHS:
***********************************************************************
NN = 1000; tMax = 50;
x0 = Table[RandomReal[{0, 1}], {i, 1, NN}];
(*x0=Table[1,{i,1,NN}];*)
AMatr = Table[RandomReal[{0, 1}], {i, 1, NN}, {j, 1, NN}];
Equations = x'[t] == - x[t]/(1 + 800 (AMatr.x[t])^2/NN^2);
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}]
{2.672, Null}
***************************************************************************
If I had taken this model from the very beginning, I would have no
problems! But I also appreciate the elegant solution by Oliver with
f[x_List]=...
Here is the non-vectorized version of this program, something I was
using before on many occasions:
**************************************************************************
NN = 100; 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}];
AMatr = Table[RandomReal[{0, 1}], {i, 1, NN}, {j, 1, NN}];
Equations = Table[x[i]'[t] == -x[i][t]/(1 + 800 Sum[AMatr[[i, j]] 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}}, PlotRange -> {0, 1}]
Out[7]= {20.921, Null}
**************************************************************************
Note that here NN=100 since NN=1000 crashes on my laptop because of
the lack of memory. The execution time 21 for the non-vectorized
version with NN=100 is longer than the execution time 2.7 for the
vectorized version with N=1000.
The bottom line is that vectorization of systems of ODEs solved by
NDSolve brings a great advantage both in speed and in memory usage.
Thank you for an excellent support again, it was extremely important
for me!
Dmitry