MathGroup Archive 2006

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

Search the Archive

Re: ndsolve & compile

  • To: mathgroup at smc.vnet.net
  • Subject: [mg64977] Re: ndsolve & compile
  • From: rknapp at wolfram.com
  • Date: Fri, 10 Mar 2006 05:15:08 -0500 (EST)
  • References: <dulsk2$397$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Alas, but you have run into the dark side of  giving vector equations.
Currently, Mathematica cannot deduce the Jacobian or its sparse
structure from the vector equations.

When you add diffusion to our equations, they become stiff, so
Mathematica switches
to implicit methods, which require a Jacobian evaluation.  If the
structure of the Jacobian is not known, this is very expensive.

The built in PDE solver, when it does the spatial discretization, keeps
track of the structure of the Jacobian.  This makes a huge difference
in the solution time:

In[65]:=
Timing[NDSolve[{D[n[t, x], t] ==
     r*n[t, x]*(1 - n[t, x]/k) - c*n[t, x]*
       (p[t, x]/(n[t, x] + kp)) + d*D[n[t, x], x, x],
    D[p[t, x], t] ==
     c*n[t, x]*(p[t, x]/(n[t, x] + kp)) +
      d*D[p[t, x], x, x], n[0, x] ==
     UnitStep[x - 0.4] - UnitStep[x - 0.6], n[t,0] == n[t,1],
    p[0, x] == 0.3, p[t,0] == p[t,1]}, {n, p}, {t, 0, tmax}, {x, 0, 1},

   Method -> {"MethodOfLines",
     "SpatialDiscretization" -> {"TensorProductGrid",
       "MinPoints" -> xmax, "MaxPoints" -> xmax,
       "DifferenceOrder" -> 2}}]]

>From In[65]:=
NDSolve::eerri:
   Warning: Estimated initial error on the specified
     spatial grid in the direction of independent
     variable x exceeds prescribed error tolerance.
     More...

Out[65]=
{0.36 Second, {{n ->

     InterpolatingFunction[{{0., 100.},

       {..., 0., 1., ...}}, <>],

    p -> InterpolatingFunction[{{0., 100.},

       {..., 0., 1., ...}}, <>]}}}

This is exactly the same system as you specified by the vector system
and is about as fast on my machine as the vector system with d = 0.
Note that the message is issued because of the discontinuous initial
condition.


There is a way to specity the structure of the Jacobian through the
Jacobian option to NDSolve, but I don't advise doing this unless you
combine n and p into a single vector  so there can be no possible
confusion about variable ordering.

Another alternative that will be avialable in the future is to use
stabilized Runge Kutta methods that are currently under development.
These are explicit methods that can extend their region of stability
for equations just like these.


Rob Knapp
Wolfram Research


  • Prev by Date: Re: Fastest method for comparing overlapping times in random time series
  • Next by Date: Re: Mathematica and Education
  • Previous by thread: Re: ndsolve & compile
  • Next by thread: Distinguished logarithm, branch cuts, etc.