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