ndsolve & compile
- To: mathgroup at smc.vnet.net
- Subject: [mg64885] ndsolve & compile
- From: Christopher Klausmeier <klausme1 at msu.edu>
- Date: Mon, 6 Mar 2006 05:01:52 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
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 ("[mg63909]
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