Re: ndsolve & compile
- To: mathgroup at smc.vnet.net
- Subject: [mg64915] Re: ndsolve & compile
- From: rknapp at wolfram.com
- Date: Tue, 7 Mar 2006 06:12:21 -0500 (EST)
- References: <duh3bs$5ut$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
A big part of the speedup I obtained from Alberto Verga's code was by
using Mathematica's listable arithmetic.
However, there is another potential pitfall with precompiling a
function that you have
fallen into, so I'll show the advantages both of working around this,
and of using listable
arithmetic.
I couldn't find a tmax in your messge for timing comparison, so I just
took the setting of
it to be tmax = 100. With your parameter settings, your compiled
function took 5.8 Seconds on my machine versus 0.312 seconds to handle
the final case you show
where you NDSolve process the equations (and automatically compile the
function)
The problem with the CompiledFunction cf you have generated is that it
was compiled
without knowledge of the values of the parameters. To resolve this, it
has to call the
Mathematica evaluator each time a parameter is used to get the current
value -- without
this, the CompiledFunction would not be correct if you changed values
of the parameters.
To work around this, you want to compile the function with the specific
values of the
parameters set. This is a little tricky, since Compile is HoldAll, so
parameters do not get evaluated. A convenient way to work around this
is to use the way Mathematica's
pattern matcher substitutes values to get specific values of the
parameters into the function.
Here is a simple example.
In[1]:=
a=1;
c1=Compile[{{n,_Integer}},Module[{h=N[Pi]/n},Table[Sin[a h
i],{i,0,n}]]]
Out[2]=
CompiledFunction[{n},
N[Pi]
Module[{h = -----}, Table[Sin[a h i], {i, 0, n}]],
n
-CompiledCode-]
Note that a is represented as a symbol in the output of the
CompiledFunction.
Its value is checked each time through the loop, which is very
expensive
In[3]:=
Timing[c1[10^6];]
Out[3]=
{5.656 Second, Null}
Here, I've jsut moved the evaluation of the parameter outside the loop
-- and the result is
the same but much faster.
In[4]:=
a=1;
c2=Compile[{{n,_Integer}},Module[{h=a N[Pi]/n},Table[Sin[h
i],{i,0,n}]]]
Out[5]=
CompiledFunction[{n},
N[Pi]
Module[{h = a -----}, Table[Sin[h i], {i, 0, n}]],
n
-CompiledCode-]
In[6]:=
Timing[c2[10^6];]
Out[6]=
{0.235 Second, Null}
Setting up something to make sure evaluation of parameters is outside
of the loop
can be cumbersome, and for your DEs would be inconvenient.
When you define a function like f[x_] := x^2, Mathematica's pattern
matcher looks for
a match of f[x_] (f[whatever] will be a match), then substitutes the
values of named
patterns (x) in the right hand side *before* evaluating. This is
independent of evaluation
so the substitution is done whether or not there are HoldAll functions
in the right
hand side
In[7]:=
CompileSin[a_] := Compile[{{n,_Integer}},Module[{h=
N[Pi]/n},Table[Sin[a h i],{i,0,n}]]];
c3 = CompileSin[a]
Out[8]=
CompiledFunction[{n},
N[Pi]
Module[{h$ = -----},
n
Table[Sin[1 h$ i], {i, 0, n}]], -CompiledCode-]
When I call CompileSin[a] with a set to 1, the rhs becomes
Compile[{{n,_Integer}},Module[{h= N[Pi]/n},Table[Sin[1 h i],{i,0,n}]]
evaluating to the result you see above. This is fast since it does not
have to check the value of a.
In[9]:=
Timing[c3[10^6];]
Out[9]=
{0.234 Second, Null}
When you use substitution to set parameters, you have to keep in mind
that
whenever you change parameter values, you have to recompile:
In[10]:=
n = 2000;hn = 1001;
In[11]:=
{c1[n][[hn]], c2[n][[hn]], c3[n][[hn]]}
Out[11]=
{1., 1., 1.}
In[12]:=
a = 2;
{c1[n][[hn]], c2[n][[hn]], c3[n][[hn]]}
Out[13]=
-16 -16
{1.22461 10 , 1.22461 10 , 1.}
c3 was compiled for a specific value of a, so it does not pick up the
change in a.
Here is how you can use this method to compile your function.
In[11]:=
CompileRHS[r_, k_, kp_, c_, m_] := Compile[{{n,_Real,1},{p,_Real,1}},
Module[{dndt,dpdt,nx=Length[n]},Transpose[Table[dndt=r*n[[x]]*(1-n[[
x]]/k)-c*n[[x]]*p[[x]]/(n[[x]]+kp);
dpdt=c*n[[x]]*p[[x]]/(n[[x]]+kp)-m*p[[x]];
{dndt,dpdt},{x,1,nx}]]]];
In[12]:=
cf1 = CompileRHS[r, k, kp, c, m];
g1[{n_?VectorQ,p_?VectorQ}]:=cf1[n,p];
First[Timing[sol=First[NDSolve[{np'[t]\[Equal]g1[np[
t]],np[0]\[Equal]{Table[0.69807,{x,
1,nx}],Table[0.32324,{x,1,nx}]}},np,{t,0,tmax}]]]]
Out[14]=
0.312 Second
which is about as fast as just giving NDSolve the DEs -- this is to be
expected since
NDSolve automatically compiles the equations in a form a bit more
efficient than using Part.
In general, using Part inside a loop is not the fastest way to evaluate
expressions
in Mathematica. Mathematica's listable arithmetic is highly optimized
and is
usually faster.
Here I compile the function to use listable arithmetic:
In[15]:=
CompileRHSListable[r_, k_, kp_, c_, m_] :=
Compile[{{n,_Real,1},{p,_Real,1}},
{
r*n*(1 - n/k) - c*n*p/(n + kp),
c*n*p/(n + kp) - m*p
}];
cfL = CompileRHSListable[r, k, kp, c, m];
gL[{n_?VectorQ,p_?VectorQ}]:=cfL[n,p];
First[Timing[sol=
First[NDSolve[{np'[t]\[Equal]gL[np[
t]],np[0]\[Equal]{Table[0.69807,{x,1,
nx}],Table[0.32324,{x,1,nx}]}},np,{t,0,tmax}]]]]
Out[18]=
0.188 Second
Getting an extra factor of 2 out of the situtation.
Note that the advantages of using the compiler diminish when you are
using listable
arithmetic since the arithmetic operations are the same whether you
compile or not.
You can give the equations in a vector form
In[19]:=
First[Timing[sol = NDSolve[
{{n'[t] \[Equal] r*n[t]*(1 - n[t]/k) - c*n[t]*p[t]/(n[t] + kp),
p'[t] \[Equal] c*n[t]*p[t]/(n[t] + kp) - m*p[t]},
{n[0] \[Equal] Table[0.69807,{x,1,nx}],
p[0] \[Equal] Table[0.32324,{x,1,nx}]}},{n,p},
{t,0,tmax}]]]
Out[19]=
0.157 Second
This is a little faster since the separation of n and p does not have
to be handled
by the pattern matcher.
Note that in this example, the _?VectorQ is not needed since with
symbolic n[t] and p[t]
the rhs of the equations do not evaluate.
Sorry for the long reply, but this issues are not trivial and require
some patience to understand.
Rob Knapp
Wolfram Research
Christopher Klausmeier wrote:
> Hi folks --
>
> Recently Rob Knapp helped greatly speed up Alberto Verga's code to
> solve a large set of ODEs by using Compile with NDSolve ("
> Re: NDSolve useless?" of 21 Jan 2006). I was hoping to see such a
> speedup myself, but got a big slowdown instead. Any ideas what is
> different?
>
> thanks -- Chris
>
> In[213]:=
> r=1;
> k=2;
> kp=1;
> c=1;
> m=0.1;
>
> nx=128;
>
> cf=Compile[{{n,_Real,1},{p,_Real,1}},
> Module[{dndt,dpdt,nx=Length[n]},Transpose[
> Table[
> dndt=r*n[[x]]*(1-n[[x]]/k)-c*n[[x]]*p[[x]]/(n[[x]]+kp);
> dpdt=c*n[[x]]*p[[x]]/(n[[x]]+kp)-m*p[[x]];
> {dndt,dpdt},{x,1,nx}]]]
> ];
>
> g[{n_?VectorQ,p_?VectorQ}]:=cf[n,p];
>
> First[Timing[
> sol=First[NDSolve[{np'[t] == g[np[t]],np[0] == {Table[0.69807,
> {x,1,nx}],
> Table[0.32324,{x,1,nx}]}},np,{t,0,tmax}]]]]
> Out[221]=
> 26.2489 Second
>
>
> In[222]:=
> eqns=Flatten[Table[{
>
> n[x]'[t]\[Equal]
> r*n[x][t]*(1-n[x][t]/k)-c*n[x][t]*p[x][t]/(n[x][t]+kp),
> p[x]'[t]\[Equal]c*n[x][t]*p[x][t]/(n[x][t]+kp)-m*p[x][t]},
> {x,1,
> nx}]];
>
> ics=Flatten[Table[{n[x][0]\[Equal]0.69807,p[x][0]\[Equal]0.32324},{x,
> 1,nx}]];
>
> unks=Flatten[Table[{n[x],p[x]},{x,1,nx}]];
>
> First[Timing[NDSolve[Flatten[Join[eqns,ics]],unks,{t,0,tmax}]]]
> Out[225]=
> 1.50484 Second
> --
> Christopher Klausmeier
> Kellogg Biological Station
> Michigan State University
> Hickory Corners MI 49060
>
> Web: http://www.kbs.msu.edu/Faculty/Klausmeier/Index.htm
> Email: klausme1 at msu.edu
> Phone: (269) 671-4987