Vector Runge-Kutta ODE solver with compilation?
DmitryG
Inspired by the new capability of Mathematica 8 to compile pieces of
user code with its own or external C compiler, I am trying to
implement the simple 4th-order Runge-Kutta solver with fixed step in a
vector form, that is, suitable for solving an arbitrary number of
equations.
The solver without compilation, the idea of which is taken from
http://www.jasoncantarella.com, is working fine:
(* Definition *)
RK4[f_,x0_,t0_,tMax_,n_] := Module[{h,K1,K2, K3,
K4,x=x0,SolList={{t0,x0}}},
h = (tMax - t0)/n;
Do[
t=t0+k h;
K1 = h f[t,x];
K2 = h f[t+(1/2) h, x +(1/2) K1];
K3 = h f[t+(1/2) h, x + (1/2) K2];
K4 = h f[t+h, x + K3];
x=x + (1/6) K1 +(1/3) K2 +(1/3) K3 +(1/6) K4;
SolList=Append[SolList,{t,x}]
,{k,1,n}];
SolList
];
(* Testing for a system of two equations *)
F[t_,x_] := {x[[2]],1 - x[[1]] + 4/(x[[1]]^3)};
t0=AbsoluteTime[];
Solution=RK4[ F,{1.0,1.0},0.0,200.0,5000];
AbsoluteTime[]-t0
ListLinePlot[Take[Solution[[All,2]],100],PlotMarkers-
>Automatic,PlotStyle->{Red}]
Out[55]= 0.5468750
The output above is execution time. The original code uses NestList
instead of my Do cycle but the execution time is the same.
Unfortunately, the compiled version of this code does not work as
expected:
(* Definition *)
RK4Comp =
Compile[{{f, _Real, 1}, {x0, _Real, 1}, {t0, _Real}, {tMax, _Real},
{n, _Integer}},
Module[{h, K1, K2, K3, K4, SolList = {{t0, x0}}, x = x0, t},
h = (tMax - t0)/n;
Do[
t = t0 + k h;
K1 = h f[t, x];
K2 = h f[t + (1/2) h, x + (1/2) K1];
K3 = h f[t + (1/2) h, x + (1/2) K2];
K4 = h f[t + h, x + K3];
x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
SolList = Append[SolList, {t, x}]
, {k, 1, n}];
SolList
](*,CompilationTarget->"C"*)
]
(* Testing for a system of two ODEs *)
F[t_, x_] := {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)};
t0 = AbsoluteTime[];
Solution = RK4Comp[ F, {1.0, 1.0}, 0.0, 200.0, 5000];
AbsoluteTime[] - t0
ListLinePlot[Take[Solution[[All, 2]], 100], PlotMarkers -> Automatic,
PlotStyle -> {Red}]
During evaluation of In[57]:= CompiledFunction::cfta: Argument F at
position 1 should be a rank 1 tensor of machine-size real numbers. >>
Out[61]= 0.5312500
Mathematica complains and seemingly proceeds without compilation
because execution time is the same.
Anybody has an idea of what is happening and how it can be corected? I
believe the problem should be relevant for many.
Thank you for your time,
Dmitry
Oliver Ruebenkoenig
Dear DmityG,
This has the problem that it is not working with PackedArrays which will
not deliver the results with good performance.
Developer`PackedArrayQ@Solution
False
Also, note that you store the time but remove it with Solution[[All,2]].
>
> The output above is execution time. The original code uses NestList
> instead of my Do cycle but the execution time is the same.
>
> Unfortunately, the compiled version of this code does not work as
> expected:
>
It works as expected, but, for example, the input f needs to be a real
vector - not a function definition. There is a difference if you have a
function that returns a vector or a (compiled) function that expects a
vector as input.
That's what this message states. The fact the this still runs is because
the compiled code will always revert to the original code when the
compiled code can not evaluate.
> Out[61]= 0.5312500
>
> Mathematica complains and seemingly proceeds without compilation
> because execution time is the same.
>
> Anybody has an idea of what is happening and how it can be corected? I
> believe the problem should be relevant for many.
>
Here is how I'd do it:
Since this is Mathematica, I'd write a function that will generate me the
compiled function like so:
makeCompRK[fIn_] := With[{f = fIn},
Compile[{{x0, _Real,
1}, {t0, _Real}, {tMax, _Real}, {n, _Integer}},
Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
h = (tMax - t0)/n;
SolList = Table[{0., 0.}, {n + 1}];
SolList[[1]] = x0;
Do[
t = t0 + k h;
K1 = h f[t, x];
K2 = h f[t + (1/2) h, x + (1/2) K1];
K3 = h f[t + (1/2) h, x + (1/2) K2];
K4 = h f[t + h, x + K3];
x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
SolList[[k]] = x;
, {k, 2, n + 1}];
SolList], CompilationTarget -> "C"]];
Note, that I have thrown away the input f and instead inject the code for
the ode into the compiled function. (That's the With[] stuff). Also, I do
not store t - if you really need that (perhaps because you'd like to
implement an adaptive RK) I can show you how to do that too. Also, I
replaced the AppendTo with a Table and set the parts of it directly.
(Small note: ConstantArray will not to the trick as it does not compile as
you could see with CompilePrint)
OK, then
f = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)}];
RK4Comp = makeCompRK[f];
t0 = AbsoluteTime[];
Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 500000];
AbsoluteTime[] - t0
runs in about 0.75 sec. on my laptop - that is for n = 500000 (and not
n=5000)
Developer`PackedArrayQ@Sol
True
Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 5000];
ListLinePlot[Take[Sol, 100], PlotMarkers -> Automatic,
PlotStyle -> {Red}]
I hope this helps,
Oliver
Oliver Ruebenkoenig
DmitryG,
one thing I forgot to mention is
tutorial/NDSolvePlugIns
that will show you how you can then use your own method and plug that into
NDSolve.
Oliver
> The output above is execution time. The original code uses NestList
> instead of my Do cycle but the execution time is the same.
>
> Unfortunately, the compiled version of this code does not work as
> expected:
>
> Out[61]= 0.5312500
>
> Mathematica complains and seemingly proceeds without compilation
> because execution time is the same.
>
> Anybody has an idea of what is happening and how it can be corected? I
> believe the problem should be relevant for many.
>
DmitryG
s
> you could see with CompilePrint)
>
> OK, then
>
> f = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)}];
> RK4Comp = makeCompRK[f];
>
> t0 = AbsoluteTime[];
> Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 500000];
> AbsoluteTime[] - t0
>
> runs in about 0.75 sec. on my laptop - that is for n = 500000 (and not
> n=5000)
>
> Developer`PackedArrayQ@Sol
> True
>
> Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 5000];
> ListLinePlot[Take[Sol, 100], PlotMarkers -> Automatic,
> PlotStyle -> {Red}]
>
> I hope this helps,
>
> Oliver
>
>
>
>
>
>
>
> > Thank you for your time,
>
> > Dmitry
Dear Oliver,
thank you for your help! Your code works on my computer.
I see there is a whole bunch of problems and I had absolutely no
chance to solve any of them myself! Mathematica documentation in Help
and online gives only simple examples of compilations where
compilation is probably not necessary. Vector RK4 is an important
generic real-life problem and, I think, its implementation has to be
included in Mathematica documentation.
I see the idea to compile RK4 with the particular system of ODEs. I
was trying to compile a RK4 that takes any system of ODEs but it did
not work. I believe, in any case, RK4 compiled with a particular
system of ODEs should be faster.
Of course, I do need t in the solution. My attempt to return to my
method of building the solution list with SolList = Append[SolList,
{t, x}] does not work. It seems compiler does not understand lists of
a more complicated structure. Also modification of your method of
initiatingSolList to zeros and then defining SolList[[k]]={t, x} does
not work. In fact, whatever I do to put t into SolList, it does not
compile. Could you, help with this, please?
Best,
Dmitry
DmitryG
Oliver,
This looks nice and works but is very slow. What is the step size h in
the fixed-step RK4 algorithm? Here is my working version for the same
equations:
**************************************************************************
CRK4[]["Step"[rhs_, t_, h_, y_, yp_]] := Module[{k0, k1, k2, k3},
k0 = h yp;
k1 = h rhs[t + h/2, y + k0/2];
k2 = h rhs[t + h/2, y + k1/2];
k3 = h rhs[t + h, y + k2];
{h, (k0 + 2 k1 + 2 k2 + k3)/6}]
CRK4[___]["DifferenceOrder"] := 4
CRK4[___]["StepMode"] := Fixed (* What is "Fixed"?? Mathematica does
not know this. What is the step h?? *)
Equations = {x1'[t] == x2[t], x2'[t] == 1 - x1[t] + 4/(x1[t]^3)};
IniConds = {x1[0] == 1, x2[0] == 1};
tMax = 200;
tt0 = AbsoluteTime[];
fixed = NDSolve[Join[Equations, IniConds], {x1, x2}, {t, 0, tMax},
Method -> CRK4, MaxSteps -> 1000000]
AbsoluteTime[] - tt0
x1t[t_] := x1[t] /. fixed[[1]];
x2t[t_] := x2[t] /. fixed[[1]];
ParametricPlot[{x1t[t], x2t[t]}, {t, 0, tMax}]
Out[61]=3.0000000
***************************************************************************=
*
Here, for tMax=500, the execution time is 3.00 and the plot shows that
the precision is insufficient because of a too large h. For larger
tMax it becomes dramatic.
Best,
Dmitry
DmitryG
> On 21 Feb., 06:08, Oliver Ruebenkoenig <ruebe...@wolfram.com> wrote:
>
>
>
>
>
>
>
>
>
> > Dear DmityG,
>
> > On Mon, 21 Feb 2011, DmitryG wrote:
> > > Dear All,
>
> > > Inspired by the new capability of Mathematica 8 to compile pieces of
> > > user code with its own or external C compiler, I am trying to
l
> > not deliver the results with good performance.
>
> > Developer`PackedArrayQ@Solution
> > False
>
.
>
> > > The output above is execution time. The original code uses NestList
> > > instead of my Do cycle but the execution time is the same.
>
> > > Unfortunately, the compiled version of this code does not work as
> > > expected:
>
> > It works as expected, but, for example, the input f needs to be a real
> > vector - not a function definition. There is a difference if you have a
> > function that returns a vector or a (compiled) function that expects a
> > vector as input.
>
> > > (* Definition *)
> > > RK4Comp =
,
> > > {n, _Integer}},
> > > Module[{h, K1, K2, K3, K4, SolList = {{t0, x0}}, x = x0, t},
>
> > > h = (tMax - t0)/n;
>
> > > Do[
> > > t = t0 + k h;
> > > K1 = h f[t, x];
> > > K2 = h f[t + (1/2) h, x + (1/2) K1];
> > > K3 = h f[t + (1/2) h, x + (1/2) K2];
> > > K4 = h f[t + h, x + K3];
> > > x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
>
> > > SolList = Append[SolList, {t, x}]
> > > , {k, 1, n}];
> > > SolList
> > > ](*,CompilationTarget->"C"*)
> > > ]
>
> > > (* Testing for a system of two ODEs *)
> > > F[t_, x_] := {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)};
> > > t0 = AbsoluteTime[];
> > > Solution = RK4Comp[ F, {1.0, 1.0}, 0.0, 200.0, 5000];
> > > AbsoluteTime[] - t0
> > > ListLinePlot[Take[Solution[[All, 2]], 100], PlotMarkers -> Automatic,
> > > PlotStyle -> {Red}]
>
> > > During evaluation of In[57]:= CompiledFunction::cfta: Argument F at
> > > position 1 should be a rank 1 tensor of machine-size real numbers. >>
>
e
> > the compiled code will always revert to the original code when the
> > compiled code can not evaluate.
>
> > > Out[61]= 0.5312500
>
> > > Mathematica complains and seemingly proceeds without compilation
> > > because execution time is the same.
>
I
> > > believe the problem should be relevant for many.
>
> > Here is how I'd do it:
>
he
> > compiled function like so:
>
> > makeCompRK[fIn_] := With[{f = fIn},
> > Compile[{{x0, _Real,
> > 1}, {t0, _Real}, {tMax, _Real}, {n, _Integer}},
> > Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
> > h = (tMax - t0)/n;
> > SolList = Table[{0., 0.}, {n + 1}];
> > SolList[[1]] = x0;
> > Do[
> > t = t0 + k h;
> > K1 = h f[t, x];
> > K2 = h f[t + (1/2) h, x + (1/2) K1];
> > K3 = h f[t + (1/2) h, x + (1/2) K2];
> > K4 = h f[t + h, x + K3];
> > x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
> > SolList[[k]] = x;
> > , {k, 2, n + 1}];
> > SolList], CompilationTarget -> "C"]];
>
or
> > the ode into the compiled function. (That's the With[] stuff). Also, I =
do
> > not store t - if you really need that (perhaps because you'd like to
> > implement an adaptive RK) I can show you how to do that too. Also, I
> > replaced the AppendTo with a Table and set the parts of it directly.
a=
> s
> > you could see with CompilePrint)
>
> > OK, then
>
> > f = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)}];
> > RK4Comp = makeCompRK[f];
>
> > t0 = AbsoluteTime[];
> > Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 500000];
> > AbsoluteTime[] - t0
>
Well, after many tries I have written a compilable vector RK4 ODE
solver that stored time t in the output list:
***********************************************************************
(* Definition *)
makeCompRK[fIn_] :=
With[{f = fIn},
Compile[{{x0, _Real, 1}, {t0, _Real}, {tMax, _Real}, {n, _Integer}},
Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
h = (tMax - t0)/n;
SolList = Table[Prepend[x0, t0], {n + 1}];
Do[t = t0 + k h;
K1 = h f[t, x];
K2 = h f[t + (1/2) h, x + (1/2) K1];
K3 = h f[t + (1/2) h, x + (1/2) K2];
K4 = h f[t + h, x + K3];
x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
SolList[[k + 1]] = Prepend[x, t ]; (* This line did not work for
some reason, be careful! *)
, {k, 1, n}];
SolList
], CompilationTarget -> "C"]]
(* Calculation *)
F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)}];
RK4Comp = makeCompRK[F];
x0 = {1.0, 1.0}; t0 = 0; tMax = 200; n = 5000;
tt0 = AbsoluteTime[];
Sol = RK4Comp[x0, t0, tMax, n];
AbsoluteTime[] - tt0
Developer`PackedArrayQ@Sol
(* Plotting *)
nTake = 200;
ListLinePlot[Take[Sol[[All, {2, 3}]], nTake],
PlotMarkers -> Automatic, PlotStyle -> {Red}]
ListLinePlot[{Take[Sol[[All, {1, 2}]], nTake],
Take[Sol[[All, {1, 3}]], nTake]}, PlotMarkers -> Automatic,
PlotStyle -> {Blue, Green}]
***************************************************************************
My idea was to extend the output list SolList as {x[[1]],x[[2]],...} -
> {t,x[[1]],x[[2],...}. I had a lot of problems with the command
SolList[[k + 1]] = Prepend[x, t ]
that gave errors and non-compilation for no apparent reason but then
suddenly worked. This format of data storage is convenient to extract
individual curves x[[i]][t] as plotting shows.
However, storing t leads to doubling of the execution time in
comparison to the method without saving t. This means, probably, that
SolList[[k + 1]] = Prepend[x, t ]
consumes so much time as everything else! So I would recommend not to
store the time that can be generated after the problem has been
solved.
Best,
Dmitry
P.S.: A problem remains for me why the same method used as a plugin to
NDSolve is so much slower and less acurate.
Oliver Ruebenkoenig
> On 21 Feb., 06:08, Oliver Ruebenkoenig <ruebe...@wolfram.com> wrote:
>> DmitryG,
>>
>> one thing I forgot to mention is
>>
>> tutorial/NDSolvePlugIns
>>
>> that will show you how you can then use your own method and plug that into
>> NDSolve.
>>
>> Oliver
>>
>>
>>
>>
>>
>>
>>
>> On Mon, 21 Feb 2011, DmitryG wrote:
>>> Dear All,
>>
>>> Inspired by the new capability of Mathematica 8 to compile pieces of
>>> user code with its own or external C compiler, I am trying to
>>> implement the simple 4th-order Runge-Kutta solver with fixed step in a
>>> The output above is execution time. The original code uses NestList
>>> instead of my Do cycle but the execution time is the same.
>>
>>> Unfortunately, the compiled version of this code does not work as
>>> expected:
>>
>>> (* Definition *)
>>> RK4Comp =
>>> Compile[{{f, _Real, 1}, {x0, _Real, 1}, {t0, _Real}, {tMax, _Real},
>>> {n, _Integer}},
>>> Module[{h, K1, K2, K3, K4, SolList = {{t0, x0}}, x = x0, t},
>>
>>> h = (tMax - t0)/n;
>>
>>> Do[
>>> t = t0 + k h;
>>> K1 = h f[t, x];
>>> K2 = h f[t + (1/2) h, x + (1/2) K1];
>>> K3 = h f[t + (1/2) h, x + (1/2) K2];
>>> K4 = h f[t + h, x + K3];
>>> x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
>>
>>> SolList = Append[SolList, {t, x}]
>>> , {k, 1, n}];
>>> SolList
>>> ](*,CompilationTarget->"C"*)
>>> ]
>>
>>> (* Testing for a system of two ODEs *)
>>> F[t_, x_] := {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)};
>>> t0 = AbsoluteTime[];
>>> Solution = RK4Comp[ F, {1.0, 1.0}, 0.0, 200.0, 5000];
>>> AbsoluteTime[] - t0
>>> ListLinePlot[Take[Solution[[All, 2]], 100], PlotMarkers -> Automatic,
>>> PlotStyle -> {Red}]
>>
>>> During evaluation of In[57]:= CompiledFunction::cfta: Argument F at
>>> position 1 should be a rank 1 tensor of machine-size real numbers. >>
>>
>>> Out[61]= 0.5312500
>>
>>> Mathematica complains and seemingly proceeds without compilation
>>> because execution time is the same.
>>
>>> Anybody has an idea of what is happening and how it can be corected? I
>>> believe the problem should be relevant for many.
>>
>>> Thank you for your time,
>>
>>> Dmitry
>
Dmitry,
> Oliver,
>
> This looks nice and works but is very slow. What is the step size h in
> the fixed-step RK4 algorithm? Here is my working version for the same
> equations:
>
> **************************************************************************
>
> CRK4[]["Step"[rhs_, t_, h_, y_, yp_]] := Module[{k0, k1, k2, k3},
> k0 = h yp;
> k1 = h rhs[t + h/2, y + k0/2];
> k2 = h rhs[t + h/2, y + k1/2];
> k3 = h rhs[t + h, y + k2];
> {h, (k0 + 2 k1 + 2 k2 + k3)/6}]
> CRK4[___]["DifferenceOrder"] := 4
> CRK4[___]["StepMode"] := Fixed (* What is "Fixed"?? Mathematica does
> not know this. What is the step h?? *)
>
> Equations = {x1'[t] == x2[t], x2'[t] == 1 - x1[t] + 4/(x1[t]^3)};
> IniConds = {x1[0] == 1, x2[0] == 1};
> tMax = 200;
>
> tt0 = AbsoluteTime[];
> fixed = NDSolve[Join[Equations, IniConds], {x1, x2}, {t, 0, tMax},
> Method -> CRK4, MaxSteps -> 1000000]
> AbsoluteTime[] - tt0
You could use StartingStepSize -> 0.04 to control the step size.
>
> x1t[t_] := x1[t] /. fixed[[1]];
> x2t[t_] := x2[t] /. fixed[[1]];
>
> ParametricPlot[{x1t[t], x2t[t]}, {t, 0, tMax}]
>
> Out[61]=3.0000000
>
> ***************************************************************************=
> *
>
Just to mention this:
AbsoluteTiming[
var = NDSolve[Join[Equations, IniConds], {x1, x2}, {t, 0, tMax},
Method -> "ExplicitRungeKutta"]
]
is quite efficient and uses adaptive step sizes.
Oliver
Oliver Ruebenkoenig
> On 21 Feb., 06:08, Oliver Ruebenkoenig <ruebe...@wolfram.com> wrote:
>> Dear DmityG,
>> This has the problem that it is not working with PackedArrays which will
>> not deliver the results with good performance.
>>
>> Developer`PackedArrayQ@Solution
>> False
>>
>> Also, note that you store the time but remove it with Solution[[All,2]].
>>
>>
>>
>>> The output above is execution time. The original code uses NestList
>>> instead of my Do cycle but the execution time is the same.
>>
>>> Unfortunately, the compiled version of this code does not work as
>>> expected:
>>
>> It works as expected, but, for example, the input f needs to be a real
>> vector - not a function definition. There is a difference if you have a
>> function that returns a vector or a (compiled) function that expects a
>> vector as input.
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> That's what this message states. The fact the this still runs is because
>> the compiled code will always revert to the original code when the
>> compiled code can not evaluate.
>>
>>> Out[61]= 0.5312500
>>
>>> Mathematica complains and seemingly proceeds without compilation
>>> because execution time is the same.
>>
>>> Anybody has an idea of what is happening and how it can be corected? I
>>> believe the problem should be relevant for many.
>>
>> Here is how I'd do it:
>>
>> Since this is Mathematica, I'd write a function that will generate me the
>> compiled function like so:
>>
>> makeCompRK[fIn_] := With[{f = fIn},
>> Compile[{{x0, _Real,
>> 1}, {t0, _Real}, {tMax, _Real}, {n, _Integer}},
>> Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
>> h = (tMax - t0)/n;
>> SolList = Table[{0., 0.}, {n + 1}];
>> SolList[[1]] = x0;
>> Do[
>> t = t0 + k h;
>> K1 = h f[t, x];
>> K2 = h f[t + (1/2) h, x + (1/2) K1];
>> K3 = h f[t + (1/2) h, x + (1/2) K2];
>> K4 = h f[t + h, x + K3];
>> x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
>> SolList[[k]] = x;
>> , {k, 2, n + 1}];
>> SolList], CompilationTarget -> "C"]];
>>
>> Note, that I have thrown away the input f and instead inject the code for
>> the ode into the compiled function. (That's the With[] stuff). Also, I do
>> not store t - if you really need that (perhaps because you'd like to
>> implement an adaptive RK) I can show you how to do that too. Also, I
>> replaced the AppendTo with a Table and set the parts of it directly.
>> (Small note: ConstantArray will not to the trick as it does not compile a=
> s
>> you could see with CompilePrint)
>>
>> OK, then
>>
>> f = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)}];
>> RK4Comp = makeCompRK[f];
>>
>> t0 = AbsoluteTime[];
>> Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 500000];
>> AbsoluteTime[] - t0
>>
>> runs in about 0.75 sec. on my laptop - that is for n = 500000 (and not
>> n=5000)
>>
>> Developer`PackedArrayQ@Sol
>> True
>>
>> Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 5000];
>> ListLinePlot[Take[Sol, 100], PlotMarkers -> Automatic,
>> PlotStyle -> {Red}]
>>
>> I hope this helps,
>>
>> Oliver
>>
>>
>>
>>
>>
>>
>>
>>> Thank you for your time,
>>
>>> Dmitry
>
> Dear Oliver,
>
> thank you for your help! Your code works on my computer.
Hello Dmitry,
>
> I see there is a whole bunch of problems and I had absolutely no
> chance to solve any of them myself! Mathematica documentation in Help
> and online gives only simple examples of compilations where
> compilation is probably not necessary. Vector RK4 is an important
> generic real-life problem and, I think, its implementation has to be
> included in Mathematica documentation.
>
> I see the idea to compile RK4 with the particular system of ODEs. I
> was trying to compile a RK4 that takes any system of ODEs but it did
> not work. I believe, in any case, RK4 compiled with a particular
> system of ODEs should be faster.
>
> Of course, I do need t in the solution. My attempt to return to my
> method of building the solution list with SolList = Append[SolList,
> {t, x}] does not work. It seems compiler does not understand lists of
> a more complicated structure. Also modification of your method of
> initiatingSolList to zeros and then defining SolList[[k]]={t, x} does
> not work. In fact, whatever I do to put t into SolList, it does not
> compile. Could you, help with this, please?
The compiler returns tensors. So you have two choices. One is to return a
list of this structure
{{t0,x0,y0},{t1,x1,y1},....} which should be the first choice.
Sometimes, however, you have to revert to another mechanism: Say that t is
not real but integer then instead of
cf = Compile[{{x, _Real, 0}}, {{1, 2. + x, 3.}, {2, 3. + x, 4.}}]
cf[1.]
you could do something like
ll
cf2 = Compile[{{x, _Real, 0}},
ll = {1, 2}; {{2. + x, 3.}, {3. + x, 4.}}];
cf2[1.]
{{3., 3.}, {4., 4.}}
ll
{1,2}
Hope this helps,
Oliver
>
> Best,
>
> Dmitry
>
>
DmitryG
> On 21 Feb., 19:29, DmitryG <einsch...@gmail.com> wrote:
>
>
>
>
>
>
>
>
>
> > On 21 Feb., 06:08, Oliver Ruebenkoenig <ruebe...@wolfram.com> wrote:
>
f=
> or
> > > the ode into the compiled function. (That's the With[] stuff). Also, =
I =
> do
> > > not store t - if you really need that (perhaps because you'd like to
> > > implement an adaptive RK) I can show you how to do that too. Also, I
> > > replaced the AppendTo with a Table and set the parts of it directly.
le=
> a=
> > s
> > > you could see with CompilePrint)
>
> > > OK, then
>
> > > f = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)}];
> > > RK4Comp = makeCompRK[f];
>
> > > t0 = AbsoluteTime[];
> > > Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 500000];
> > > AbsoluteTime[] - t0
>
no=
> t
> > > n=5000)
>
> > > Developer`PackedArrayQ@Sol
> > > True
>
> > > Sol = RK4Comp[{1.0, 1.0}, 0.0, 200.0, 5000];
> > > ListLinePlot[Take[Sol, 100], PlotMarkers -> Automatic,
> > > PlotStyle -> {Red}]
>
> > > I hope this helps,
>
> > > Oliver
>
> > > > Thank you for your time,
>
> > > > Dmitry
>
>
> > thank you for your help! Your code works on my computer.
>
> > chance to solve any of them myself! Mathematica documentation in Help
> > and online gives only simple examples of compilations where
> > compilation is probably not necessary. Vector RK4 is an important
> > generic real-life problem and, I think, its implementation has to be
> > included in Mathematica documentation.
>
> > I see the idea to compile RK4 with the particular system of ODEs. I
> > was trying to compile a RK4 that takes any system of ODEs but it did
> > not work. I believe, in any case, RK4 compiled with a particular
> > system of ODEs should be faster.
>
> > Of course, I do need t in the solution. My attempt to return to my
> > method of building the solution list with SolList = Append[SolList,
> > {t, x}] does not work. It seems compiler does not understand lists of
> > a more complicated structure. Also modification of your method of
> > initiatingSolList to zeros and then defining SolList[[k]]={t, x} does
> > not work. In fact, whatever I do to put t into SolList, it does not
> > compile. Could you, help with this, please?
>
>
> > Dmitry
>
> Well, after many tries I have written a compilable vector RK4 ODE
> solver that stored time t in the output list:
>
> ***********************************************************************
> (* Definition *)
> With[{f = fIn},
> Compile[{{x0, _Real, 1}, {t0, _Real}, {tMax, _Real}, {n, _Integer}},
> Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
> h = (tMax - t0)/n;
> K1 = h f[t, x];
> K2 = h f[t + (1/2) h, x + (1/2) K1];
> K3 = h f[t + (1/2) h, x + (1/2) K2];
> K4 = h f[t + h, x + K3];
> x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
ork for
> some reason, be careful! *)
> SolList
>
> (* Calculation *)
> F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/(x[[1]]^3)}];
> RK4Comp = makeCompRK[F];
> x0 = {1.0, 1.0}; t0 = 0; tMax = 200; n = 5000;
> tt0 = AbsoluteTime[];
> Sol = RK4Comp[x0, t0, tMax, n];
> AbsoluteTime[] - tt0
> Developer`PackedArrayQ@Sol
>
> (* Plotting *)
> nTake = 200;
> PlotMarkers -> Automatic, PlotStyle -> {Red}]
>
> PlotStyle -> {Blue, Green}]
>
> *************************************************************************=
**
>
> My idea was to extend the output list SolList as {x[[1]],x[[2]],...} -
>
> > {t,x[[1]],x[[2],...}. I had a lot of problems with the command
>
> SolList[[k + 1]] = Prepend[x, t ]
>
> that gave errors and non-compilation for no apparent reason but then
> suddenly worked. This format of data storage is convenient to extract
> individual curves x[[i]][t] as plotting shows.
>
> However, storing t leads to doubling of the execution time in
> comparison to the method without saving t. This means, probably, that
>
> SolList[[k + 1]] = Prepend[x, t ]
>
> consumes so much time as everything else! So I would recommend not to
> store the time that can be generated after the problem has been
> solved.
>
> Best,
>
> Dmitry
>
> P.S.: A problem remains for me why the same method used as a plugin to
> NDSolve is so much slower and less acurate.
Dear Oliver and All interested,
Today I tried to apply the above very efficient vector RK4 method with
compilation to my real problems that I was solving with NDSolve
before. Unfortunately, it has proven to be more complicated than it
seemed. To illustrate the difficulties, I will use the same working
model as above and will try to generalize it. So, the working code is
*************************************************
(* Definition of the RK4 solver *)
makeCompRK[fIn_] :=
With[{f = fIn},
Compile[{{x0, _Real, 1}, {t0, _Real}, {tMax, _Real}, {n,
_Integer}},
Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
h = (tMax - t0)/n;
SolList = Table[x0, {n + 1}];
Do[t = t0 + k h;
K1 = h f[t, x];
K2 = h f[t + (1/2) h, x + (1/2) K1];
K3 = h f[t + (1/2) h, x + (1/2) K2];
K4 = h f[t + h, x + K3];
x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
SolList[[k + 1]] = x
, {k, 1, n}];
SolList], CompilationTarget -> "C"]];
(* Calculation *)
F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
RK4Comp = makeCompRK[F];
x0 = {1.0, 1.0}; t0 = 0; tMax = 200; n = 5000;
tt0 = AbsoluteTime[];
Sol = RK4Comp[x0, t0, tMax, n];
AbsoluteTime[] - tt0
Developer`PackedArrayQ@Sol
(* Plotting *)
nTake = 200;
ListLinePlot[Take[Sol, nTake], PlotMarkers -> Automatic, PlotStyle ->
{Red}] (* Trajectory *)
(* Time dependences *)
x1List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 1]]}, {k, 0, n}];
x2List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 2]]}, {k, 0, n}];
ListLinePlot[{Take[x1List, nTake], Take[x2List, nTake]},
PlotMarkers -> Automatic, PlotStyle -> {Blue, Green}]
***************************************************************************=
**************************
In this example, the vector of ODE's RHS' is defined by
F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
In real life, this is something more complicated, a result of a many-
step calculation. One can try to generalize the code as
RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3};
F = Function[{t, x}, RHS];
(that should be equivalent to the above) and then substitute RHS by
whatever you have for your problem. However, Mathematica complains and
produces crap results. My suspicion is that it is only F =
Function[{t, x}, RHS]; that is injected into the compiled RK4 code
whereas RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; is not injected.
However, I cannot be sure because it was never explained.
Further idea is to define F as a procedure and put everything inside
it. A toy example is
F = Function[{t, x}, RHS = {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; RHS];
This does not compile, although then retreating to the non-compilation
mode produces correct results. My experiments show that you can only
use a single expression in the body of F that does not use any
definitions. For instance,
F = Function[{t, x}, Table[-x[[i]]^i, {i, 1, 2}]];
compiles and works fine but this is, again, no more than a toy
example.
To summarize, I am very excited by the prospect to speed up my
computations by a factor about 100 but I don't see how to generalize
the approach for real-life problems.
Oliver Ruebenkoenig
[...]
Dear Dmitry,
>
> Dear Oliver,
>
> thank you for your help! Your code works on my computer.
>
> I see there is a whole bunch of problems and I had absolutely no
> chance to solve any of them myself! Mathematica documentation in Help
> and online gives only simple examples of compilations where
> compilation is probably not necessary. Vector RK4 is an important
> generic real-life problem and, I think, its implementation has to be
> included in Mathematica documentation.
>
you find more examples on how to extend NDSolve with real-life methods
here:
LibraryLink/tutorial/Numerical
Oliver
[...]
Oliver Ruebenkoenig
On Wed, 23 Feb 2011, DmitryG wrote:
[....]
Dmitry,
your friend is - in princial -
F = Function[{t, x}, Evaluate[RHS]]
That has a small problem that we used Part in the RHS and Part gives a
message that you try to access elements that the _symbol_ x does not have.
One way to work around this is
RHS = Quiet[{x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
F = Function[{t, x}, Evaluate[RHS]];
If you evaluate F you'll see that the RHS in now in the function body. If
you look at the F you posted, there the body is RHS, which is not what
you want. All this hast to do with
Attributes[Function]
HoldAll.
Oliver
DmitryG
;
> > x2List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 2]]}, {k, 0, n}]=
;
> > ListLinePlot[{Take[x1List, nTake], Take[x2List, nTake]},
> > PlotMarkers -> Automatic, PlotStyle -> {Blue, Green}]
Thank you, Oliver, this was a great help.
Now, I am trying to introduce sampling of the results to reduce the
number of output points. This is also very important. However, there
is a new problem. Here is the code:
******************************************************************************
(* Vector RK4 routine *)
makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
{ttMax,_Real},{n,_Integer},{nsamp,_Integer}},
Module[{h,K1,K2,K3,K4,SolList,x=x0,t,ksamp},
h=(ttMax-t0)/n;
SolList=Table[x0,{nsamp+1}];
Do[t=t0+k h; ksamp=k nsamp/n;
K1=h f[t,x];
K2=h f[t+(1/2) h,x+(1/2) K1];
K3=h f[t+(1/2) h,x+(1/2) K2];
K4=h f[t+h,x+K3];
x=x+(1/6) K1+(1/3) K2+(1/3) K3+(1/6) K4;
If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
,{k,1,n}];
SolList],CompilationTarget->"C"]];
(* Calculation *)
RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]; (* Quiet to suppress
complaints *)
F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
RHS *)
RK4Comp=makeCompRK[F];
x0={1.0,1.0}; t0=0; tMax; NSteps=5000; NSampling0;
tt0=AbsoluteTime[];
Sol=RK4Comp[x0,t0,tMax,NSteps,NSampling];
AbsoluteTime[]-tt0
Developer`PackedArrayQ@Sol
******************************************************************************
The problem is in the line
If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
In the first execution, Mathematica complains
Compile::cpintlt: "ksamp$+1 at position 2 of SolList$[[ksamp$+1]]
should be either a nonzero integer or a vector of nonzero integers;
evaluation will use the uncompiled function."
I do not understand this. Still Mathematica compiles and the results
are correct.
In the second execution, Mathematica complains more and does not
compile, although the results are still correct.
If I use the command
If[IntegerQ[kSampling],SolList[[kSampling+1]]=x
Mathematica complains again:
CompiledFunction::cfse: "Compiled expression False should be a machine-
size real number)."
does not compile, but the results are correct.
What is wrong here? I believe I am doing everything in a natural way.
But compiling is too tricky..
Dmitry
Leonid Shifrin
Just a side remark: if using delayed assignment is acceptable for RHS,
another way to avoid
error messages (and generally the evaluation of RHS) would be
In[48]:=
RHS:={x[[2]],1-x[[1]]+4/x[[1]]^3};
Function@@Prepend[Hold[RHS]/.OwnValues[RHS],Unevaluated[{t,x}]]
Out[49]= Function[{t,x},{x[[2]],1-x[[1]]+4/x[[1]]^3}]
I found this idiom useful on several occasions.
Cheers,
Leonid
On Wed, Feb 23, 2011 at 2:25 PM, Oliver Ruebenkoenig
<rueb...@wolfram.com>wrote:
> > x1List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 1]]}, {k, 0, n}];
> > x2List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 2]]}, {k, 0, n}];
> > ListLinePlot[{Take[x1List, nTake], Take[x2List, nTake]},
> > PlotMarkers -> Automatic, PlotStyle -> {Blue, Green}]
> >
> ***************************************************************************=
DmitryG
];
>
> > > This does not compile, although then retreating to the non-compilatio=
n
> > > mode produces correct results. My experiments show that you can only
> > > use a single expression in the body of F that does not use any
> > > definitions. For instance,
>
> > > F = Function[{t, x}, Table[-x[[i]]^i, {i, 1, 2}]];
>
> > > compiles and works fine but this is, again, no more than a toy
> > > example.
>
> > > To summarize, I am very excited by the prospect to speed up my
> > > computations by a factor about 100 but I don't see how to generalize
> > > the approach for real-life problems.
>
> Thank you, Oliver, this was a great help.
>
> Now, I am trying to introduce sampling of the results to reduce the
> number of output points. This is also very important. However, there
> is a new problem. Here is the code:
>
** ***
> (* Vector RK4 routine *)
> makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
> {ttMax,_Real},{n,_Integer},{nsamp,_Integer}},
> Module[{h,K1,K2,K3,K4,SolList,x=x0,t,ksamp},
> h=(ttMax-t0)/n;
> SolList=Table[x0,{nsamp+1}];
> Do[t=t0+k h; ksamp=k nsamp/n;
> K1=h f[t,x];
> K2=h f[t+(1/2) h,x+(1/2) K1];
> K3=h f[t+(1/2) h,x+(1/2) K2];
> K4=h f[t+h,x+K3];
> x=x+(1/6) K1+(1/3) K2+(1/3) K3+(1/6) K4;
>
> If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>
> ,{k,1,n}];
> SolList],CompilationTarget->"C"]];
>
> (* Calculation *)
> RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]; (* Quiet to suppress
> complaints *)
> F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
> RHS *)
> RK4Comp=makeCompRK[F];
;
> tt0=AbsoluteTime[];
> Sol=RK4Comp[x0,t0,tMax,NSteps,NSampling];
> AbsoluteTime[]-tt0
> Developer`PackedArrayQ@Sol
** ***
>
> The problem is in the line
>
> If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>
> In the first execution, Mathematica complains
>
> Compile::cpintlt: "ksamp$+1 at position 2 of SolList$[[ksamp$+1]]
> should be either a nonzero integer or a vector of nonzero integers;
> evaluation will use the uncompiled function."
>
> I do not understand this. Still Mathematica compiles and the results
> are correct.
>
> In the second execution, Mathematica complains more and does not
> compile, although the results are still correct.
>
> If I use the command
>
> If[IntegerQ[kSampling],SolList[[kSampling+1]]=x
>
> Mathematica complains again:
>
> CompiledFunction::cfse: "Compiled expression False should be a machine-
> size real number)."
>
> does not compile, but the results are correct.
>
> What is wrong here? I believe I am doing everything in a natural way.
> But compiling is too tricky..
>
> Dmitry
The command that works is
If[FractionalPart[ksamp] == 0, SolList[[IntegerPart[ksamp] + 1]] x];
(* Although assignment is done for ksamp integer, IntegerPart is
needed! *)
Commands including the test
IntegerQ[ksamp]
do not work.
Dmitry
Oliver Ruebenkoenig
>>> F = Function[{t, x}, RHS = {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; RHS];
>>
>>> This does not compile, although then retreating to the non-compilation
>>> mode produces correct results. My experiments show that you can only
>>> use a single expression in the body of F that does not use any
>>> definitions. For instance,
>>
>>> F = Function[{t, x}, Table[-x[[i]]^i, {i, 1, 2}]];
>>
>>> compiles and works fine but this is, again, no more than a toy
>>> example.
>>
>>> To summarize, I am very excited by the prospect to speed up my
>>> computations by a factor about 100 but I don't see how to generalize
>>> the approach for real-life problems.
>
> Thank you, Oliver, this was a great help.
>
> Now, I am trying to introduce sampling of the results to reduce the
> number of output points. This is also very important. However, there
> is a new problem. Here is the code:
>
> ******************************************************************************
> (* Vector RK4 routine *)
> makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
> {ttMax,_Real},{n,_Integer},{nsamp,_Integer}},
> Module[{h,K1,K2,K3,K4,SolList,x=x0,t,ksamp},
> h=(ttMax-t0)/n;
> SolList=Table[x0,{nsamp+1}];
> Do[t=t0+k h; ksamp=k nsamp/n;
> K1=h f[t,x];
> K2=h f[t+(1/2) h,x+(1/2) K1];
> K3=h f[t+(1/2) h,x+(1/2) K2];
> K4=h f[t+h,x+K3];
> x=x+(1/6) K1+(1/3) K2+(1/3) K3+(1/6) K4;
>
> If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>
> ,{k,1,n}];
> SolList],CompilationTarget->"C"]];
>
> (* Calculation *)
> RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]; (* Quiet to suppress
> complaints *)
> F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
> RHS *)
> RK4Comp=makeCompRK[F];
> x0={1.0,1.0}; t0=0; tMax; NSteps=5000; NSampling0;
> tt0=AbsoluteTime[];
> Sol=RK4Comp[x0,t0,tMax,NSteps,NSampling];
> AbsoluteTime[]-tt0
> Developer`PackedArrayQ@Sol
> ******************************************************************************
>
> The problem is in the line
>
> If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>
> In the first execution, Mathematica complains
>
> Compile::cpintlt: "ksamp$+1 at position 2 of SolList$[[ksamp$+1]]
> should be either a nonzero integer or a vector of nonzero integers;
> evaluation will use the uncompiled function."
Of what type is ksamp?
ksamp = k nsamp/n
For example
Head[5 3/2]
In worst case this is a rational and you check that with
FractionalPart. But Part needs an integer at position 2. With using Floor
you can tell the compiler that ksamp is always going to be an integer.
SolList[[Floor[ksamp] + 1]]
Oliver
DmitryG
}]=
> > ;
> >>> x2List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 2]]}, {k, 0, n=
}]=
> > ;
> >>> ListLinePlot[{Take[x1List, nTake], Take[x2List, nTake]},
> >>> PlotMarkers -> Automatic, PlotStyle -> {Blue, Green}]
****** > > **************************
>
> >>> In this example, the vector of ODE's RHS' is defined by
>
> >>> F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
>
> >>> In real life, this is something more complicated, a result of a many-
> >>> step calculation. One can try to generalize the code as
>
> >>> RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3};
> >>> F = Function[{t, x}, RHS];
>
> >> Dmitry,
>
> >> your friend is - in princial -
>
> >> F = Function[{t, x}, Evaluate[RHS]]
>
> >> That has a small problem that we used Part in the RHS and Part gives a
ave.
>
> >> One way to work around this is
>
> >> RHS = Quiet[{x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
> >> F = Function[{t, x}, Evaluate[RHS]];
>
If
> >> you look at the F you posted, there the body is RHS, which is not what
> >> you want. All this hast to do with
>
> >> Attributes[Function]
>
> >> HoldAll.
>
> >> Oliver
>
> >>> (that should be equivalent to the above) and then substitute RHS by
d
> >>> produces crap results. My suspicion is that it is only F =
> >>> Function[{t, x}, RHS]; that is injected into the compiled RK4 code
> >>> whereas RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; is not injected.
> >>> However, I cannot be sure because it was never explained.
>
> >>> Further idea is to define F as a procedure and put everything inside
> >>> it. A toy example is
>
];
>
> >>> This does not compile, although then retreating to the non-compilatio=
n
> >>> mode produces correct results. My experiments show that you can only
> >>> use a single expression in the body of F that does not use any
> >>> definitions. For instance,
>
> >>> F = Function[{t, x}, Table[-x[[i]]^i, {i, 1, 2}]];
>
> >>> compiles and works fine but this is, again, no more than a toy
> >>> example.
>
> >>> To summarize, I am very excited by the prospect to speed up my
> >>> computations by a factor about 100 but I don't see how to generalize
> >>> the approach for real-life problems.
>
> > Thank you, Oliver, this was a great help.
>
> > Now, I am trying to introduce sampling of the results to reduce the
> > number of output points. This is also very important. However, there
> > is a new problem. Here is the code:
>
**** ***
> > (* Vector RK4 routine *)
> > makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
> > {ttMax,_Real},{n,_Integer},{nsamp,_Integer}},
> > Module[{h,K1,K2,K3,K4,SolList,x=x0,t,ksamp},
> > h=(ttMax-t0)/n;
> > SolList=Table[x0,{nsamp+1}];
> > Do[t=t0+k h; ksamp=k nsamp/n;
> > K1=h f[t,x];
> > K2=h f[t+(1/2) h,x+(1/2) K1];
> > K3=h f[t+(1/2) h,x+(1/2) K2];
> > K4=h f[t+h,x+K3];
> > x=x+(1/6) K1+(1/3) K2+(1/3) K3+(1/6) K4;
>
> > If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>
> > ,{k,1,n}];
> > SolList],CompilationTarget->"C"]];
>
> > (* Calculation *)
> > RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]; (* Quiet to suppress
> > complaints *)
> > F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
> > RHS *)
> > RK4Comp=makeCompRK[F];
g0;
> > tt0=AbsoluteTime[];
> > Sol=RK4Comp[x0,t0,tMax,NSteps,NSampling];
> > AbsoluteTime[]-tt0
> > Developer`PackedArrayQ@Sol
**** ***
Thank you for the response, Oliver! I have already fugured it out.
There is now a problem with compilation in C in the case of a system
of N equations. For N above 10, it takes a long time to compile, and
the compilation time is dramatically increasing with N, so that one
cannot solve any real-life problems with large N using compilation in
C. On the other hand, using Mathematca's own compiler is non-
problematic. In both cases the solution is correct. Here is the code:
***************************************************************************=
***************
(* Definition *)
makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
{tMax,_Real},{n,_Integer}},Module[{h,K1,K2,K3,K4,SolList,x=x0,t},
h=(tMax-t0)/n;
SolList=Table[x0,{n+1}];
Do[t=t0+k h;
K1=h f[t,x];
K2=h f[t+(1/2) h,x+(1/2) K1];
K3=h f[t+(1/2) h,x+(1/2) K2];
K4=h f[t+h,x+K3];
x=x+(1/6) K1+(1/3) K2+(1/3) K3+(1/6) K4;
SolList[[k+1]]=x
,{k,1,n}];
SolList],CompilationTarget->"C"]];
(* Calculation *)
NN=30;
RHS=Quiet[Table[-x[[i]]/(1+Sum[x[[j]],{j,1,NN}]^2),{i,1,NN}]]; (*
Quiet to suppress complaints *)
F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
RHS *)
tt0=AbsoluteTime[];
Timing[RK4Comp=makeCompRK[F];]
AbsoluteTime[]-tt0
Print[];
x0=Table[RandomReal[{0,1}],{i,1,NN}]; t0=0; tMax0; n=500;
tt0=AbsoluteTime[];
Sol=RK4Comp[x0,t0,tMax,n];
AbsoluteTime[]-tt0
Print["Compilation: ",Developer`PackedArrayQ@Sol]
(* Plotting time dependences *)
x1List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,1]]},{k,0,n}];
x2List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,2]]},{k,0,n}];
x3List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,3]]},{k,0,n}];
ListLinePlot[{x1List,x2List,x3List},PlotMarkers->Automatic,PlotStyle-
>{Blue,Green,Red},PlotRange->{0,1}]
Compile::nogen: A library could not be generated from the compiled
function. >>
Out[4]= {11.172,Null}
Out[5]= 142.0625000
Out[10]= 0.0468750
Compilation: True
***************************************************************************=
**********************************
Whereas execution is fast, 0.0468750, compilation takes a long time. I
measure both Timing and AbsoluteTime. The former, I believe, is the
time used by Mathematica, whereas the latter is the total time used by
Mathematica and my C compiler, Microsoft Visual C++ 2010 Express. Both
times are long. Presumably, Mathematica generates some weird code that
gives a hard time to C compiler. Here I also have an error message
concerning a library but I haven' had this message as I ran this
program before.
Any help?
Dmitry
DmitryG
> On 25 Feb., 06:38, Oliver Ruebenkoenig <ruebe...@wolfram.com> wrote:
>
>
>
>
>
>
>
> > On Thu, 24 Feb 2011, DmitryG wrote:
> > > On 23 Feb., 06:25, Oliver Ruebenkoenig <ruebe...@wolfram.com> wrote:
> > >> On Wed, 23 Feb 2011, DmitryG wrote:
>
> > >> [....]
>
> > >>> Dear Oliver and All interested,
>
ith
> > >>> compilation to my real problems that I was solving with NDSolve
> > >>> before. Unfortunately, it has proven to be more complicated than it
> > >>> seemed. To illustrate the difficulties, I will use the same working
is
>
> > >>> *************************************************
>
> > >>> (* Definition of the RK4 solver *)
> > >>> makeCompRK[fIn_] :=
> > >>> With[{f = fIn},
> > >>> Compile[{{x0, _Real, 1}, {t0, _Real}, {tMax, _Real}, {n,
> > >>> _Integer}},
> > >>> Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
> > >>> h = (tMax - t0)/n;
> > >>> SolList = Table[x0, {n + 1}];
> > >>> Do[t = t0 + k h;
> > >>> K1 = h f[t, x];
> > >>> K2 = h f[t + (1/2) h, x + (1/2) K1];
> > >>> K3 = h f[t + (1/2) h, x + (1/2) K2];
> > >>> K4 = h f[t + h, x + K3];
> > >>> x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
> > >>> SolList[[k + 1]] = x
>
> > >>> , {k, 1, n}];
> > >>> SolList], CompilationTarget -> "C"]];
>
> > >>> (* Calculation *)
> > >>> F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
> > >>> RK4Comp = makeCompRK[F];
0;
> > >>> tt0 = AbsoluteTime[];
> > >>> Sol = RK4Comp[x0, t0, tMax, n];
> > >>> AbsoluteTime[] - tt0
> > >>> Developer`PackedArrayQ@Sol
>
> > >>> (* Plotting *)
> > >>> nTake = 200;
->
> > >>> {Red}] (* Trajectory *)
>
> > >>> (* Time dependences *)
n=
> }]=
> > > ;
> > >>> x2List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 2]]}, {k, 0,=
n=
> }]=
> > > ;
> > >>> ListLinePlot[{Take[x1List, nTake], Take[x2List, nTake]},
> > >>> PlotMarkers -> Automatic, PlotStyle -> {Blue, Green}]
**=
>
> ****** > > **************************
>
>
>
>
>
>
>
>
>
>
>
> > >>> In this example, the vector of ODE's RHS' is defined by
>
> > >>> F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
>
y-
> > >>> step calculation. One can try to generalize the code as
>
> > >>> RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3};
> > >>> F = Function[{t, x}, RHS];
>
> > >> Dmitry,
>
> > >> your friend is - in princial -
>
> > >> F = Function[{t, x}, Evaluate[RHS]]
>
a
> > >> message that you try to access elements that the _symbol_ x does not=
h=
> ave.
>
> > >> One way to work around this is
>
> > >> RHS = Quiet[{x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
> > >> F = Function[{t, x}, Evaluate[RHS]];
>
y.=
> If
> > >> you look at the F you posted, there the body is RHS, which is not wh=
at
> > >> you want. All this hast to do with
>
> > >> Attributes[Function]
>
> > >> HoldAll.
>
> > >> Oliver
>
> > >>> (that should be equivalent to the above) and then substitute RHS by
an=
> d
> > >>> produces crap results. My suspicion is that it is only F =
> > >>> Function[{t, x}, RHS]; that is injected into the compiled RK4 code
> > >>> whereas RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; is not injected.
> > >>> However, I cannot be sure because it was never explained.
>
e
> > >>> it. A toy example is
>
HS=
> ];
>
> > >>> This does not compile, although then retreating to the non-compilat=
io=
> n
> > >>> mode produces correct results. My experiments show that you can onl=
y
> > >>> use a single expression in the body of F that does not use any
> > >>> definitions. For instance,
>
> > >>> F = Function[{t, x}, Table[-x[[i]]^i, {i, 1, 2}]];
>
> > >>> compiles and works fine but this is, again, no more than a toy
> > >>> example.
>
> > >>> To summarize, I am very excited by the prospect to speed up my
e
> > >>> the approach for real-life problems.
>
> > > Thank you, Oliver, this was a great help.
>
> > > Now, I am trying to introduce sampling of the results to reduce the
> > > number of output points. This is also very important. However, there
> > > is a new problem. Here is the code:
>
**=
> **** ***
> > > (* Vector RK4 routine *)
> > > makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
> > > {ttMax,_Real},{n,_Integer},{nsamp,_Integer}},
> > > Module[{h,K1,K2,K3,K4,SolList,x=x0,t,ksamp},
> > > h=(ttMax-t0)/n;
> > > SolList=Table[x0,{nsamp+1}];
> > > Do[t=t0+k h; ksamp=k nsamp/n;
> > > K1=h f[t,x];
> > > K2=h f[t+(1/2) h,x+(1/2) K1];
> > > K3=h f[t+(1/2) h,x+(1/2) K2];
> > > K4=h f[t+h,x+K3];
> > > x=x+(1/6) K1+(1/3) K2+(1/3) K3+(1/6) K4;
>
> > > If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>
> > > ,{k,1,n}];
> > > SolList],CompilationTarget->"C"]];
>
> > > (* Calculation *)
> > > RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]; (* Quiet to suppress
> > > complaints *)
of
> > > RHS *)
> > > RK4Comp=makeCompRK[F];
in=
> g0;
> > > tt0=AbsoluteTime[];
> > > Sol=RK4Comp[x0,t0,tMax,NSteps,NSampling];
> > > AbsoluteTime[]-tt0
> > > Developer`PackedArrayQ@Sol
**=
> **** ***
>
> > > The problem is in the line
>
> > > If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>
> > > In the first execution, Mathematica complains
>
> > > Compile::cpintlt: "ksamp$+1 at position 2 of SolList$[[ksamp$+1]]
> > > should be either a nonzero integer or a vector of nonzero integers;
> > > evaluation will use the uncompiled function."
>
> > Of what type is ksamp?
>
> > ksamp = k nsamp/n
>
> > For example
>
> > Head[5 3/2]
>
> > In worst case this is a rational and you check that with
or
> > you can tell the compiler that ksamp is always going to be an integer.
>
> > SolList[[Floor[ksamp] + 1]]
>
> > Oliver
>
> > > I do not understand this. Still Mathematica compiles and the results
> > > are correct.
>
> > > In the second execution, Mathematica complains more and does not
> > > compile, although the results are still correct.
>
> > > If I use the command
>
> > > If[IntegerQ[kSampling],SolList[[kSampling+1]]=x
>
> > > Mathematica complains again:
>
e-
> > > size real number)."
>
> > > does not compile, but the results are correct.
>
> > > What is wrong here? I believe I am doing everything in a natural way.
> > > But compiling is too tricky..
>
> > > Dmitry
>
> Thank you for the response, Oliver! I have already fugured it out.
>
> There is now a problem with compilation in C in the case of a system
> of N equations. For N above 10, it takes a long time to compile, and
> the compilation time is dramatically increasing with N, so that one
> cannot solve any real-life problems with large N using compilation in
> C. On the other hand, using Mathematca's own compiler is non-
> problematic. In both cases the solution is correct. Here is the code:
>
** =
** =
> **********************************
>
> Whereas execution is fast, 0.0468750, compilation takes a long time. I
> measure both Timing and AbsoluteTime. The former, I believe, is the
> time used by Mathematica, whereas the latter is the total time used by
> Mathematica and my C compiler, Microsoft Visual C++ 2010 Express. Both
> times are long. Presumably, Mathematica generates some weird code that
> gives a hard time to C compiler. Here I also have an error message
> concerning a library but I haven' had this message as I ran this
> program before.
>
> Any help?
>
> Dmitry
Another major problem: If in the code above one removes
CompilationTarget->"C", it is working well until NN becomes larger
then ~1000. For larger NN, say NN00, compilation breaks because of
lack of memory on my laptop. One can see how fast the memory is
getting consumed during compilation. This seems to be a deficiency of
the current version of the Mathematica compiler that probably makes it
inappropriate for solving larger problems where compilation is really
needed. And there seems to be an even bigger problem with using
external C compiler.
Dmitry
Oliver Ruebenkoenig
> On 25 Feb., 06:38, Oliver Ruebenkoenig <ruebe...@wolfram.com> wrote:
>> On Thu, 24 Feb 2011, DmitryG wrote:
>>> On 23 Feb., 06:25, Oliver Ruebenkoenig <ruebe...@wolfram.com> wrote:
>>>> On Wed, 23 Feb 2011, DmitryG wrote:
>>
>>>> [....]
>>
>>>>> Dear Oliver and All interested,
>>
>>>>> Today I tried to apply the above very efficient vector RK4 method with
>>>>> compilation to my real problems that I was solving with NDSolve
>>>>> before. Unfortunately, it has proven to be more complicated than it
>>>>> seemed. To illustrate the difficulties, I will use the same working
>>>>> model as above and will try to generalize it. So, the working code is
>>
>>>>> *************************************************
>>
>>>>> (* Definition of the RK4 solver *)
>>>>> makeCompRK[fIn_] :=
>>>>> With[{f = fIn},
>>>>> Compile[{{x0, _Real, 1}, {t0, _Real}, {tMax, _Real}, {n,
>>>>> _Integer}},
>>>>> Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
>>>>> h = (tMax - t0)/n;
>>>>> SolList = Table[x0, {n + 1}];
>>>>> Do[t = t0 + k h;
>>>>> K1 = h f[t, x];
>>>>> K2 = h f[t + (1/2) h, x + (1/2) K1];
>>>>> K3 = h f[t + (1/2) h, x + (1/2) K2];
>>>>> K4 = h f[t + h, x + K3];
>>>>> x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
>>>>> SolList[[k + 1]] = x
>>
>>>>> , {k, 1, n}];
>>>>> SolList], CompilationTarget -> "C"]];
>>
>>>>> (* Calculation *)
>>>>> F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
>>>>> RK4Comp = makeCompRK[F];
>>>>> x0 = {1.0, 1.0}; t0 = 0; tMax = 200; n = 5000;
>>>>> tt0 = AbsoluteTime[];
>>>>> Sol = RK4Comp[x0, t0, tMax, n];
>>>>> AbsoluteTime[] - tt0
>>>>> Developer`PackedArrayQ@Sol
>>
>>>>> (* Plotting *)
>>>>> nTake = 200;
>>>>> ListLinePlot[Take[Sol, nTake], PlotMarkers -> Automatic, PlotStyle ->
>>>>> {Red}] (* Trajectory *)
>>
>>>>> (* Time dependences *)
>>>>> x1List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 1]]}, {k, 0, n=
> }]=
>>> ;
>>>>> x2List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 2]]}, {k, 0, n=
> }]=
>>> ;
>>>>> ListLinePlot[{Take[x1List, nTake], Take[x2List, nTake]},
>>>>> PlotMarkers -> Automatic, PlotStyle -> {Blue, Green}]
>>>>> *********************************************************************=
> ****** > > **************************
>>
>>>>> In this example, the vector of ODE's RHS' is defined by
>>
>>>>> F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
>>
>>>>> In real life, this is something more complicated, a result of a many-
>>>>> step calculation. One can try to generalize the code as
>>
>>>>> RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3};
>>>>> F = Function[{t, x}, RHS];
>>
>>>> Dmitry,
>>
>>>> your friend is - in princial -
>>
>>>> F = Function[{t, x}, Evaluate[RHS]]
>>
>>>> That has a small problem that we used Part in the RHS and Part gives a
>>>> message that you try to access elements that the _symbol_ x does not h=
> ave.
>>
>>>> One way to work around this is
>>
>>>> RHS = Quiet[{x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
>>>> F = Function[{t, x}, Evaluate[RHS]];
>>
>>>> If you evaluate F you'll see that the RHS in now in the function body.=
> If
>>>> you look at the F you posted, there the body is RHS, which is not what
>>>> you want. All this hast to do with
>>
>>>> Attributes[Function]
>>
>>>> HoldAll.
>>
>>>> Oliver
>>
>>>>> (that should be equivalent to the above) and then substitute RHS by
>>>>> whatever you have for your problem. However, Mathematica complains an=
> d
>>>>> produces crap results. My suspicion is that it is only F =
>>>>> Function[{t, x}, RHS]; that is injected into the compiled RK4 code
>>>>> whereas RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; is not injected.
>>>>> However, I cannot be sure because it was never explained.
>>
>>>>> Further idea is to define F as a procedure and put everything inside
>>>>> it. A toy example is
>>
>>>>> F = Function[{t, x}, RHS = {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; RHS=
> ];
>>
>>>>> This does not compile, although then retreating to the non-compilatio=
> n
>>>>> mode produces correct results. My experiments show that you can only
>>>>> use a single expression in the body of F that does not use any
>>>>> definitions. For instance,
>>
>>>>> F = Function[{t, x}, Table[-x[[i]]^i, {i, 1, 2}]];
>>
>>>>> compiles and works fine but this is, again, no more than a toy
>>>>> example.
>>
>>>>> To summarize, I am very excited by the prospect to speed up my
>>>>> computations by a factor about 100 but I don't see how to generalize
>>>>> the approach for real-life problems.
>>
>>> Thank you, Oliver, this was a great help.
>>
>>> Now, I am trying to introduce sampling of the results to reduce the
>>> number of output points. This is also very important. However, there
>>> is a new problem. Here is the code:
>>
>>> ***********************************************************************=
> **** ***
>>> (* Vector RK4 routine *)
>>> makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
>>> {ttMax,_Real},{n,_Integer},{nsamp,_Integer}},
>>> Module[{h,K1,K2,K3,K4,SolList,x=x0,t,ksamp},
>>> h=(ttMax-t0)/n;
>>> SolList=Table[x0,{nsamp+1}];
>>> Do[t=t0+k h; ksamp=k nsamp/n;
>>> K1=h f[t,x];
>>> K2=h f[t+(1/2) h,x+(1/2) K1];
>>> K3=h f[t+(1/2) h,x+(1/2) K2];
>>> K4=h f[t+h,x+K3];
>>> x=x+(1/6) K1+(1/3) K2+(1/3) K3+(1/6) K4;
>>
>>> If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>>
>>> ,{k,1,n}];
>>> SolList],CompilationTarget->"C"]];
>>
>>> (* Calculation *)
>>> RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]; (* Quiet to suppress
>>> complaints *)
>>> F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
>>> RHS *)
>>> RK4Comp=makeCompRK[F];
>>> x0={1.0,1.0}; t0=0; tMax; NSteps=5000; NSamplin=
> g0;
>>> tt0=AbsoluteTime[];
>>> Sol=RK4Comp[x0,t0,tMax,NSteps,NSampling];
>>> AbsoluteTime[]-tt0
>>> Developer`PackedArrayQ@Sol
>>> ***********************************************************************=
> **** ***
>>
>>> The problem is in the line
>>
>>> If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>>
>>> In the first execution, Mathematica complains
>>
>>> Compile::cpintlt: "ksamp$+1 at position 2 of SolList$[[ksamp$+1]]
>>> should be either a nonzero integer or a vector of nonzero integers;
>>> evaluation will use the uncompiled function."
>>
>> Of what type is ksamp?
>>
>> ksamp = k nsamp/n
>>
>> For example
>>
>> Head[5 3/2]
>>
>> In worst case this is a rational and you check that with
>> FractionalPart. But Part needs an integer at position 2. With using Floor
>> you can tell the compiler that ksamp is always going to be an integer.
>>
>> SolList[[Floor[ksamp] + 1]]
>>
>> Oliver
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>> I do not understand this. Still Mathematica compiles and the results
>>> are correct.
>>
>>> In the second execution, Mathematica complains more and does not
>>> compile, although the results are still correct.
>>
>>> If I use the command
>>
>>> If[IntegerQ[kSampling],SolList[[kSampling+1]]=x
>>
>>> Mathematica complains again:
>>
>>> CompiledFunction::cfse: "Compiled expression False should be a machine-
>>> size real number)."
>>
>>> does not compile, but the results are correct.
>>
>>> What is wrong here? I believe I am doing everything in a natural way.
>>> But compiling is too tricky..
>>
>>> Dmitry
>
> Thank you for the response, Oliver! I have already fugured it out.
>
> There is now a problem with compilation in C in the case of a system
> of N equations. For N above 10, it takes a long time to compile, and
You could use CompilePrint to see that for NN=11 this generates a some 900
line so code - quite a lot of which are Part[]. Make sure you understand
what you are doing here - you in-line this sum into the code.
> the compilation time is dramatically increasing with N, so that one
> cannot solve any real-life problems with large N using compilation in
> C. On the other hand, using Mathematca's own compiler is non-
> problematic. In both cases the solution is correct. Here is the code:
>
> ***************************************************************************=
> x0=Table[RandomReal[{0,1}],{i,1,NN}]; t0=0; tMax0; n=500;
> tt0=AbsoluteTime[];
> Sol=RK4Comp[x0,t0,tMax,n];
> AbsoluteTime[]-tt0
>
> Print["Compilation: ",Developer`PackedArrayQ@Sol]
>
> (* Plotting time dependences *)
>
> x1List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,1]]},{k,0,n}];
> x2List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,2]]},{k,0,n}];
> x3List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,3]]},{k,0,n}];
> ListLinePlot[{x1List,x2List,x3List},PlotMarkers->Automatic,PlotStyle-
>> {Blue,Green,Red},PlotRange->{0,1}]
>
>
> Compile::nogen: A library could not be generated from the compiled
> function. >>
> Out[4]= {11.172,Null}
> Out[5]= 142.0625000
>
>
> Out[10]= 0.0468750
> Compilation: True
>
> ***************************************************************************=
> **********************************
>
> Whereas execution is fast, 0.0468750, compilation takes a long time. I
> measure both Timing and AbsoluteTime. The former, I believe, is the
> time used by Mathematica, whereas the latter is the total time used by
> Mathematica and my C compiler, Microsoft Visual C++ 2010 Express. Both
> times are long. Presumably, Mathematica generates some weird code that
Weird code? - No, not really.
> gives a hard time to C compiler. Here I also have an error message
> concerning a library but I haven' had this message as I ran this
> program before.
>
> Any help?
>
> Dmitry
>
>
You have several options
1)
F = With[{NN = NN}, Compile[{{t, _Real, 0}, {x, _Real, 1}},
Table[-x[[i]]/(1 + Sum[x[[j]], {j, 1, NN}]^2), {i, 1, NN}],
CompilationTarget -> C]];
Use a compiled function and insert that into the RK4 code. This is like
another function call.
- If you look at your code - that is does a lot of computations
repetitively, the sum is computed for every i - I do not think that any
compiler can optimize that.
2) Use vector code
F = Compile[{{t, _Real, 0}, {x, _Real, 1}}, -x/(1 + Total[x]^2),
CompilationTarget -> C];
3) To save the additional function call you can the inline this piece
code.
RHS = Quiet[(-x/(1 + Total[x]^2))]
F = Function[{t, x}, Evaluate[RHS]]
NN = 30;
x0 = Table[
RandomReal[{0, 1}], {i, 1, NN}]; t0 = 0; tMax = 2000; n = 50000;
(* n = 50000 *)
tt0 = AbsoluteTime[];
Timing[RK4Comp = makeCompRK[F];]
AbsoluteTime[] - tt0
{0.104007, Null}
0.257110
tt0 = AbsoluteTime[];
Sol = RK4Comp[x0, t0, tMax, n];
AbsoluteTime[] - tt0
1.481468
Oliver
DmitryG
> You could use CompilePrint to see that for NN=11 this generates a some 900
> line so code - quite a lot of which are Part[]. Make sure you understand
> what you are doing here - you in-line this sum into the code.
>
>
......................
>
> Weird code? - No, not really.
> ...
>
Hi Oliver,
In this short program I mimick my real-life program that solves a
system of NN ODE's with the maximal coupling: the RHS of each of NN
equations contains all NN variables. Any molecular dynamics (MD)
simulation is of this kind: Acceleration of any particle depends on
positions of all particles.
This application of Mathematica's compilation is crucial but there are
two problems that I have encountered:
1) Compilation in C is prohibitively slow above NN=50 or so. In this
case Mathematica's own compiler is still doing a good job.
2) For NN=1000 or above Mathematica's compiler consumes all RAM (my
laptop has 3 GB RAM) and quits. This is very disappointing because
generating the list of RHS does not take much memory and the RK 4
routine does not consume much memory as well. Thus it would be only a
matter of computer time to solve a problem for any NN. However,
compiler generates a blown-up code (component by component, with the
vector structure abandoned) that does not fit into memory.
Thank you for suggestion to look into the compiled code with
CompilePrint. Yes, the code is very long, and its structure is the
same for the initial problem we considered, the problem with
RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]
Also in this case the compiled code contains Part[].
Is there any way to optimize this calculation or we have an intrinsic
limitation of Mathematica here? Below is the code again:
*******************************************************************************
(* Definition *)
makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
{tMax,_Real},{n,_Integer}},Module[{h,K1,K2,K3,K4,SolList,x=x0,t},
h=(tMax-t0)/n;
SolList=Table[x0,{n+1}];
Do[t=t0+k h;
K1=h f[t,x];
K2=h f[t+(1/2) h,x+(1/2) K1];
K3=h f[t+(1/2) h,x+(1/2) K2];
K4=h f[t+h,x+K3];
x=x+(1/6) K1+(1/3) K2+(1/3) K3+(1/6) K4;
SolList[[k+1]]=x
,{k,1,n}];
SolList](*,CompilationTarget->"C"*)]];
(* Calculation *)
NN=15;
Timing[RHS=Quiet[Table[-x[[i]]/(1+300Sum[x[[j]],{j,1,NN}]^2/NN^2),{i,
1,NN}]];] (* Quiet to suppress complaints *)
F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
RHS *)
Print[]
tt0=AbsoluteTime[];
Timing[RK4Comp=makeCompRK[F];]
AbsoluteTime[]-tt0
Print[];
x0=Table[RandomReal[{0,1}],{i,1,NN}]; t0=0; tMax=50; n=500;
tt0=AbsoluteTime[];
Sol=RK4Comp[x0,t0,tMax,n];
AbsoluteTime[]-tt0
Print["Compilation: ",Developer`PackedArrayQ@Sol]
(* Plotting time dependences *)
x1List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,1]]},{k,0,n}];
x2List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,2]]},{k,0,n}];
x3List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,3]]},{k,0,n}];
ListLinePlot[{x1List,x2List,x3List},PlotMarkers->Automatic,PlotStyle-
>{Blue,Green,Red},PlotRange->{0,1}]
*********************************
Thank you for your help,
Dmitry
Daniel Lichtblau
----- Original Message -----
> [...big snip...]
> > ***************************************************************************
> > ***************
> > (* Definition *)
> > makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
> > {tMax,_Real},{n,_Integer}},Module[{h,K1,K2,K3,K4,SolList,x=x0,t},
> > h=(tMax-t0)/n;
> > SolList=Table[x0,{n+1}];
> > Do[t=t0+k h;
> > K1=h f[t,x];
> > K2=h f[t+(1/2) h,x+(1/2) K1];
> > K3=h f[t+(1/2) h,x+(1/2) K2];
> > K4=h f[t+h,x+K3];
> > x=x+(1/6) K1+(1/3) K2+(1/3) K3+(1/6) K4;
> > SolList[[k+1]]=x
> >
> > ,{k,1,n}];
> > SolList],CompilationTarget->"C"]];
> >
> > (* Calculation *)
> > NN=30;
> > RHS=Quiet[Table[-x[[i]]/(1+Sum[x[[j]],{j,1,NN}]^2),{i,1,NN}]]; (*
> > Quiet to suppress complaints *)
> > F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
> > RHS *)
> >
> > tt0=AbsoluteTime[];
> > Timing[RK4Comp=makeCompRK[F];]
> > AbsoluteTime[]-tt0
> > Print[];
> >
> > x0=Table[RandomReal[{0,1}],{i,1,NN}]; t0=0; tMax0; n=
> =500;
> > tt0=AbsoluteTime[];
> > Sol=RK4Comp[x0,t0,tMax,n];
> > AbsoluteTime[]-tt0
> >
> > Print["Compilation: ",Developer`PackedArrayQ@Sol]
> >
> > (* Plotting time dependences *)
> >
> > x1List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,1]]},{k,0,n}];
> > x2List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,2]]},{k,0,n}];
> > x3List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,3]]},{k,0,n}];
> > ListLinePlot[{x1List,x2List,x3List},PlotMarkers->Automatic,PlotStyle-
> >
> > >{Blue,Green,Red},PlotRange->{0,1}]
> >
> > Compile::nogen: A library could not be generated from the compiled
> > function. >>
> > Out[4]= {11.172,Null}
> > Out[5]= 142.0625000
> >
> > Out[10]= 0.0468750
> > Compilation: True
> >
> > *************************************************************************=
> ** =
> > **********************************
> >
> > Whereas execution is fast, 0.0468750, compilation takes a long time.
> > I
> > measure both Timing and AbsoluteTime. The former, I believe, is the
> > time used by Mathematica, whereas the latter is the total time used
> > by
> > Mathematica and my C compiler, Microsoft Visual C++ 2010 Express.
> > Both
> > times are long. Presumably, Mathematica generates some weird code
> > that
> > gives a hard time to C compiler. Here I also have an error message
> > concerning a library but I haven' had this message as I ran this
> > program before.
> >
> > Any help?
> >
> > Dmitry
>
> Another major problem: If in the code above one removes
> CompilationTarget->"C", it is working well until NN becomes larger
> then ~1000. For larger NN, say NN00, compilation breaks because of
> lack of memory on my laptop. One can see how fast the memory is
> getting consumed during compilation. This seems to be a deficiency of
> the current version of the Mathematica compiler that probably makes it
> inappropriate for solving larger problems where compilation is really
> needed. And there seems to be an even bigger problem with using
> external C compiler.
>
> Dmitry
This is an unusual claim to make. The implication is something along the lines "Surely my usage is flawless, so all problems must be the fault of the underlying software."
I code a fair amount in both Mathematica and C, and that sort of claim is one I would rarely make with either.
Daniel Lichtblau
Wolfram Research
DmitryG
>
> > > Any help?
>
> > > Dmitry
>
> > Another major problem: If in the code above one removes
> > CompilationTarget->"C", it is working well until NN becomes larger
> > then ~1000. For larger NN, say NN00, compilation breaks because of
> > lack of memory on my laptop. One can see how fast the memory is
> > getting consumed during compilation. This seems to be a deficiency of
> > the current version of the Mathematica compiler that probably makes it
> > inappropriate for solving larger problems where compilation is really
> > needed. And there seems to be an even bigger problem with using
> > external C compiler.
>
> > Dmitry
>
> This is an unusual claim to make. The implication is something along the lines "Surely my usage is flawless, so all problems must be the fault of the underlying software."
>
> I code a fair amount in both Mathematica and C, and that sort of claim is one I would rarely make with either.
>
> Daniel Lichtblau
> Wolfram Research
Hi Daniel,
I am posting a working code and you can quickly run it with and
without compilation in C, to see the problem.
I do not claim that my code is flawless. It is straightforward but
maybe you could suggest how to eliminate the bottleneck in compilation
in C. I would be very grateful if you do. Currently I cannot use this
option for my problems.
Compilation with Mathematica's own compiler works fine up to NN=1200
on my standard laptop and I am considering using other computers with
more memory. BTW, Mathematica 7 on my (older) workstation also has
Compile implemented but it is working much slower than Mathematica 8
on my laptop. This is to say that I see the progress and I hope that
compilation will be improved in future releases of Mathematica.
Best,
Dmitry
Daniel Lichtblau
> makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
> {tMax,_Real},{n,_Integer}},Module[{h,K1,K2,K3,K4,SolList,x=x0,t},
> h=(tMax-t0)/n;
> SolList=Table[x0,{n+1}];
> Do[t=t0+k h;
> K1=h f[t,x];
> K2=h f[t+(1/2) h,x+(1/2) K1];
> K3=h f[t+(1/2) h,x+(1/2) K2];
> K4=h f[t+h,x+K3];
> x=x+(1/6) K1+(1/3) K2+(1/3) K3+(1/6) K4;
> SolList[[k+1]]=x
>
> ,{k,1,n}];
>
> (* Calculation *)
> NN=15;
> Timing[RHS=Quiet[Table[-x[[i]]/(1+300Sum[x[[j]],{j,1,NN}]^2/NN^2),{i,
> F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
> RHS *)
>
> tt0=AbsoluteTime[];
> Timing[RK4Comp=makeCompRK[F];]
> AbsoluteTime[]-tt0
> Print[];
>
> tt0=AbsoluteTime[];
> Sol=RK4Comp[x0,t0,tMax,n];
> AbsoluteTime[]-tt0
>
> Print["Compilation: ",Developer`PackedArrayQ@Sol]
>
> (* Plotting time dependences *)
>
> x1List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,1]]},{k,0,n}];
> x2List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,2]]},{k,0,n}];
> x3List=Table[{1.t0+(tMax-t0)k/n,Sol[[k+1,3]]},{k,0,n}];
> ListLinePlot[{x1List,x2List,x3List},PlotMarkers->Automatic,PlotStyle-
>> {Blue,Green,Red},PlotRange->{0,1}]
>
>
>
> Thank you for your help,
>
> Dmitry
>
I'll show code below that seems to be better behaved. I believe it is
functionally the same as yours but you'll need to check that. Itt is
based on some things Oliver Ruebenkoenig suggested (possibly you had not
yet seen that post at the time you wrote this). The main change is to
vectorize that RHS. I also use a Table instead of loop to fill in the
solution, but that probably is not too important.
makeCompRK[fIn_] :=
With[{f = fIn},
Compile[{{x0, _Real,
1}, {t0, _Real}, {tMax, _Real}, {n, _Integer}},
Module[{h, K1, K2, K3, K4, x = x0, t},
h = (tMax - t0)/n;
Prepend[Table[t = t0 + k h;
K1 = h f[t, x];
K2 = h f[t + (1/2.) h, x + (1/2.) K1];
K3 = h f[t + (1/2.) h, x + (1/2.) K2];
K4 = h f[t + h, x + K3];
x = x + (1/6.) K1 + (1/3.) K2 + (1/3.) K3 + (1/6.) K4,
{k, 1, n}], x0]
](*,CompilationTarget->"C"*)]];
(*Calculation*)
NN = 150;
Timing[RHS =
Quiet[-x/(1 +
300 Total[x]^2/NN^2)];] (*Quiet to suppress complaints*)
F =
Function[{t, x}, Evaluate[RHS]];(*Evaluate to use actual form of RHS*)
Out[136]= {0., Null}
tt0 = AbsoluteTime[];
Timing[RK4Comp = makeCompRK[F];]
AbsoluteTime[] - tt0
x0 =
Table[RandomReal[{0, 1}], {i, 1, NN}]; t0 = 0; tMax = 50; n = 500;
tt0 = AbsoluteTime[];
Sol = RK4Comp[x0, t0, tMax, n];
AbsoluteTime[] - tt0
Out[144]= 0.004684
Print["Compilation: ", Developer`PackedArrayQ@Sol]
(*Plotting time dependences*)
x1List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 1]]}, {k, 0,
n}];
x2List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 2]]}, {k, 0,
n}];
x3List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 3]]}, {k, 0,
n}];
ListLinePlot[{x1List, x2List, x3List}, PlotMarkers -> Automatic,
PlotStyle -> {Blue, Green, Red}, PlotRange -> {0, 1}]
With this version I had no difficulty handling NN=1500 either using C or
the virtual machine interpreter of Compile.
Daniel Lichtblau
Wolfram Research
Oliver Ruebenkoenig
> On 27 Feb., 04:43, Oliver Ruebenkoenig <ruebe...@wolfram.com> wrote:
>
>> You could use CompilePrint to see that for NN=11 this generates a some 900
>> line so code - quite a lot of which are Part[]. Make sure you understand
>> what you are doing here - you in-line this sum into the code.
>>
>>
> ......................
>>
>> Weird code? - No, not really.
>> ...
>>
Hi Dmitry,
>
> Hi Oliver,
>
> In this short program I mimick my real-life program that solves a
Why do you not send the example that you classify as real-life?
> system of NN ODE's with the maximal coupling: the RHS of each of NN
> equations contains all NN variables. Any molecular dynamics (MD)
> simulation is of this kind: Acceleration of any particle depends on
> positions of all particles.
>
One method I have seen for MD programs is to use a cut-off range that
mimics for example the Lennard-Jones potential - which reduces the
complexity of the alg.
> This application of Mathematica's compilation is crucial but there are
> two problems that I have encountered:
>
> 1) Compilation in C is prohibitively slow above NN=50 or so. In this
> case Mathematica's own compiler is still doing a good job.
>
I think the in-lining of the sum is not what you want to do here. You
should either use vectorized code or another compiled function, as shown
in a previous post.
> 2) For NN=1000 or above Mathematica's compiler consumes all RAM (my
> laptop has 3 GB RAM) and quits. This is very disappointing because
> generating the list of RHS does not take much memory and the RK 4
> routine does not consume much memory as well. Thus it would be only a
> matter of computer time to solve a problem for any NN. However,
> compiler generates a blown-up code (component by component, with the
> vector structure abandoned) that does not fit into memory.
>
No it is not that compiler that generated this - you told it to in-line that
into the code and that's what happened. - A compiler is optimized to
optimize hand written code - it can not do logical simplifications.
Would you copy and past your sum function in a C file? No, you'd write a
function and call it. In Mathematica you additionally have the option to
use vectorized code. What is wrong with a vector approach?
> Thank you for suggestion to look into the compiled code with
> CompilePrint. Yes, the code is very long, and its structure is the
> same for the initial problem we considered, the problem with
>
> RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]
>
Just the slight difference that in the initial problem NN=2.
> Also in this case the compiled code contains Part[].
>
> Is there any way to optimize this calculation or we have an intrinsic
> limitation of Mathematica here? Below is the code again:
>
The limitation is in the code you inline. It is not that you inline, but
what you inline.
Here is the code for NN=1000, tMax0 and n = 50000
makeCompRK[fIn_] :=
With[{f = fIn},
Compile[{{x0, _Real,
1}, {t0, _Real}, {tMax, _Real}, {n, _Integer}},
Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
h = (tMax - t0)/n;
SolList = Table[x0, {n + 1}];
Do[t = t0 + k h;
K1 = h f[t, x];
K2 = h f[t + (1/2) h, x + (1/2) K1];
K3 = h f[t + (1/2) h, x + (1/2) K2];
K4 = h f[t + h, x + K3];
x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
SolList[[k + 1]] = x, {k, 1, n}];
SolList], CompilationTarget -> "C"]];
NN = 1000;
RHS = Quiet[-x/(1 + 300*Total[x]^2/NN^2)]
F = Function[{t, x}, Evaluate[RHS]]
tt0 = AbsoluteTime[];
Timing[RK4Comp = makeCompRK[F];]
AbsoluteTime[] - tt0
x0 = RandomReal[{0, 1}, {NN}]; t0 = 0; tMax = 200; n = 50000;
tt0 = AbsoluteTime[];
Sol = RK4Comp[x0, t0, tMax, n];
AbsoluteTime[] - tt0
runs in about 3.7 sec on my laptop.
Oliver
DmitryG
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?
DmitryG
*********************************************************************
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[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
DrMajorBob
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[1][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[0] == x0}, x, {t, 0, tMax},
MaxSteps -> 1000000];]
{0.048119, Null}
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}]
That seems to match the following plot:
Clear[x]
initial = Table[x[i][0] == 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[1][t] /. solution[[1]];
x2t[t_] := x[2][t] /. solution[[1]];
x3t[t_] := x[3][t] /. solution[[1]];
p = Plot[{x1t[t], x2t[t], x3t[t]}, {t, 0, tMax},
PlotStyle -> {{Thick, Red}, {Thick, Green}, {Thick, Blue}}];
]
p
Bobby
Oliver Ruebenkoenig
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
Daniel Lichtblau
> 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);
>
simply evaluates to t. So the input to NDSolve is not what you want.
To retain a vector equation you might replace that with
Equations :=
x'[t] == -x[t]/(1 + 300 (ConstantArray[1, NN].x[t])^2/NN^2);
Daniel Lichtblau
Wolfram Research
DrMajorBob
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 <eins...@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
>
>
DmitryG
> 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[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}]
>
> *************************************************************************** **
> 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
I see that in my vectorized equation Total[x[t]] is misinterpreted,
possibly because Mathematica does not know x[t] is a vector.
Unfortunately I do not find in Help how to tell Mathematica x[t] is a
vector.
Dmitry
DmitryG
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
DmitryG
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
Oliver Ruebenkoenig
> On 3 Mrz., 06:00, Oliver Ruebenkoenig <ruebe...@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