MathGroup Archive 2006

[Date Index] [Thread Index] [Author Index]

Search the Archive

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





  • Prev by Date: FindRoot Secant Method
  • Next by Date: Re: mathematica to word
  • Previous by thread: FindRoot Secant Method
  • Next by thread: Re: ndsolve & compile