Stiffness - Maple Help

Numerical Solution of Stiff IVPs in Maple

Description

 • If it takes too long to solve an initial value problem (IVP) with the default numerical solver rkf45 or ck45 or with the default 'adamsfunc' option in lsode, then the IVP might be "stiff". A stiff IVP can be solved efficiently with either the default solver using the 'stiff=true' option (rosenbrock) or any of the `back-' options in lsode.
 • rosenbrock is intended for modest accuracy, so by default it uses a relative tolerance of $\frac{1}{1000}$.  You should use lsode if you need high accuracy.  By default it uses a relative tolerance of $\frac{1}{10000000}$. Numerical methods for stiff IVPs must form partial derivatives.  rosenbrock does this analytically, but if this is not possible, you can use lsode because it approximates them with finite differences.
 • Stiffness is complicated because it involves both the IVP and the method used to solve it. A stiff IVP has a solution that is very stable (some solutions that start near it converge to it very quickly) and is slowly varying (other solutions change on much faster time scales). Special methods allow stiff problems to be solved efficiently, but they can be very inefficient for problems that are not stiff.
 • Unless you have reason to believe your IVP is stiff, you should try a method for nonstiff problems first, such as the default method of dsolve[numeric] or lsode with its default option.  If that proves unsatisfactory, try the option 'stiff=true' in the call to dsolve[numeric], which invokes a Rosenbrock method. You might also try the one of the 'back-' options of lsode, which invoke backward differentiation methods (BDFs, also known as Gear's method).

Examples

A simple example of a stiff linear IVP:

 > $\mathrm{deq1}≔\mathrm{diff}\left(y\left(t\right),t\right)+10000\left(y\left(t\right)-\mathrm{cos}\left(t\right)\right)=-\mathrm{sin}\left(t\right)$
 ${\mathrm{deq1}}{≔}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){+}{10000}{}{y}{}\left({t}\right){-}{10000}{}{\mathrm{cos}}{}\left({t}\right){=}{-}{\mathrm{sin}}{}\left({t}\right)$ (1)
 > $\mathrm{ic1}≔y\left(0\right)=2$
 ${\mathrm{ic1}}{≔}{y}{}\left({0}\right){=}{2}$ (2)
 > $\mathrm{dsol1}≔\mathrm{dsolve}\left(\left\{\mathrm{deq1},\mathrm{ic1}\right\},\mathrm{numeric},\mathrm{range}=0..10,\mathrm{stiff}=\mathrm{true}\right)$
 ${\mathrm{dsol1}}{≔}{\mathbf{proc}}\left({\mathrm{x_rosenbrock}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (3)
 > $\mathrm{dsol1}\left(10.0\right)$
 $\left[{t}{=}{10.0}{,}{y}{}\left({t}\right){=}{-0.839071537030767}\right]$ (4)
 > $\mathrm{cos}\left(10.0\right)$
 ${-0.8390715291}$ (5)

Issue the following command to see a plot of the solution.

 > $\mathrm{plots}\left[\mathrm{odeplot}\right]\left(\mathrm{dsol1},\left[t,y\left(t\right)\right]\right)$

Trying to solve this IVP with the stiff=false option is too expensive. From the general solution of this equation, $\mathrm{cos}\left(t\right)+\mathrm{_C1}{ⅇ}^{-10000t}$, we see that the solution of the IVP is very stable because all solutions come together exponentially fast.  And, after an initial transient, the solution of the IVP approaches the slowly varying solution $\mathrm{cos}\left(t\right)$.

We can solve this same example with the backfull option for lsode, which by default asks for more accuracy than rosenbrock.

 > $\mathrm{dsol2}≔\mathrm{dsolve}\left(\left\{\mathrm{deq1},\mathrm{ic1}\right\},\mathrm{numeric},\mathrm{method}=\mathrm{lsode}\left[\mathrm{backfull}\right]\right)$
 ${\mathrm{dsol2}}{≔}{\mathbf{proc}}\left({\mathrm{x_lsode}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (6)
 > $\mathrm{dsol2}\left(10.0\right)$
 $\left[{t}{=}{10.0}{,}{y}{}\left({t}\right){=}{-0.839071701470174}\right]$ (7)
 > $\mathrm{cos}\left(10.0\right)$
 ${-0.8390715291}$ (8)

The van der Pol equation in relaxation oscillation is a famous example of a stiff nonlinear problem.

 > $\mathrm{\mu }≔1000$
 ${\mathrm{\mu }}{≔}{1000}$ (9)
 > $\mathrm{deq3}≔\mathrm{diff}\left(y\left(x\right),x,x\right)-\mathrm{\mu }\left(1-{y\left(x\right)}^{2}\right)\mathrm{diff}\left(y\left(x\right),x\right)+y\left(x\right)=0$
 ${\mathrm{deq3}}{≔}\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{x}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({x}\right){-}{1000}{}\left({1}{-}{{y}{}\left({x}\right)}^{{2}}\right){}\left(\frac{{ⅆ}}{{ⅆ}{x}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({x}\right)\right){+}{y}{}\left({x}\right){=}{0}$ (10)
 > $\mathrm{ic3}≔\left\{y\left(0\right)=2,\mathrm{D}\left(y\right)\left(0\right)=0\right\}$
 ${\mathrm{ic3}}{≔}\left\{{y}{}\left({0}\right){=}{2}{,}{\mathrm{D}}{}\left({y}\right){}\left({0}\right){=}{0}\right\}$ (11)
 > $\mathrm{dsol3}≔\mathrm{dsolve}\left(\left\{\mathrm{deq3}\right\}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathbf{union}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathrm{ic3},\mathrm{numeric},\mathrm{range}=0..3000,\mathrm{stiff}=\mathrm{true}\right)$
 ${\mathrm{dsol3}}{≔}{\mathbf{proc}}\left({\mathrm{x_rosenbrock}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (12)

Issue the following command to see a plot of the solution that shows regions of sharp change.

 > $\mathrm{plots}\left[\mathrm{odeplot}\right]\left(\mathrm{dsol3},\left[x,y\left(x\right)\right]\right)$

If you experiment with other initial conditions, you will see that all non-trivial solutions converge very rapidly to this limit cycle.  The IVP is stiff where the solution is slowly varying.

A stiff nonlinear system from reactor kinetics:

 > $\mathrm{deq4}≔\mathrm{diff}\left(u\left(t\right),t\right)=0.01-\left(1+\left(u\left(t\right)+1000\right)\left(u\left(t\right)+1\right)\right)\left(0.01+u\left(t\right)+v\left(t\right)\right),\mathrm{diff}\left(v\left(t\right),t\right)=0.01-\left(1+{v\left(t\right)}^{2}\right)\left(0.01+u\left(t\right)+v\left(t\right)\right)$
 ${\mathrm{deq4}}{≔}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{u}{}\left({t}\right){=}{0.01}{-}\left({1}{+}\left({u}{}\left({t}\right){+}{1000}\right){}\left({u}{}\left({t}\right){+}{1}\right)\right){}\left({0.01}{+}{u}{}\left({t}\right){+}{v}{}\left({t}\right)\right){,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{v}{}\left({t}\right){=}{0.01}{-}\left({1}{+}{{v}{}\left({t}\right)}^{{2}}\right){}\left({0.01}{+}{u}{}\left({t}\right){+}{v}{}\left({t}\right)\right)$ (13)
 > $\mathrm{ic4}≔u\left(0\right)=0,v\left(0\right)=0:$
 > $\mathrm{dsol4}≔\mathrm{dsolve}\left(\left\{\mathrm{deq4},\mathrm{ic4}\right\},\mathrm{numeric},\mathrm{method}=\mathrm{rosenbrock},\mathrm{range}=0..20\right)$
 ${\mathrm{dsol4}}{≔}{\mathbf{proc}}\left({\mathrm{x_rosenbrock}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (14)
 > $\mathrm{plots}\left[\mathrm{odeplot}\right]\left(\mathrm{dsol4},\left[t,u\left(t\right),v\left(t\right)\right]\right)$

References

 Gear, C. W. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice Hall, 1971.
 Hairer, E., and Wanner, G. Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems. Springer, 1996.
 Shampine, L. F., and Corless, R. M. "Initial Value Problems for ODEs in Problem Solving Environments." Journal of Computational and Applied Mathematics, Vol. 125 (2000): 31-40.