
Description


•

Invoking the dsolve function with the numeric option and $\mathrm{output}=\mathrm{procedurelist}$ (default), $\mathrm{output}=\mathrm{listprocedure}$ or output=operator causes dsolve to return procedure(s) to numerically solve the input ODE system.

•

Prior to returning this procedure, the input system is converted to a first order system (by introduction of new variables if necessary), then solved with respect to those first order derivatives giving a system of the form:

$x\'\left(t\right)\=\mathrm{f1}\left(t\,x\left(t\right)\,y\left(t\right)\right)$
$y\'\left(t\right)\=\mathrm{f2}\left(t\,x\left(t\right)\,y\left(t\right)\right)$
•

During the course of the numerical solution of the ODE system, it is necessary to evaluate the functions $\mathrm{f1}$, $\mathrm{f2}$ above for different values of $t$, $x$, $y$, to obtain the values of the derivatives $x\'$, $y\'$.


The option maxfun=n in the call to dsolve,numeric gives an upper limit of $n$ on the total number of evaluations performed on $\mathrm{f1}$, $\mathrm{f2}$ for any single call to the procedure returned by dsolve. Note that the evaluation of both $\mathrm{f1}$, $\mathrm{f2}$ at one set of values of $t$, $x$, $y$ is considered to be a single evaluation.

•

This option is provided as a means of limiting the amount of work performed on any single call to the solution procedure, and as a means of preventing the computation from running forever when it encounters a singularity in the solution of the ODE.

•

When the maxfun function evaluation limit is exceeded, the solution procedure halts with an error of the form Error, (in <proc>) too many function evaluations in <method>.

•

For the dverk78, lsode, and gear methods, this option is disabled by default. For the rkf45, ck45, and rosenbrock methods, the default setting is maxfun=30000, while for the classical methods, the default setting is maxfun=50000.

•

A setting of maxfun=0 disables the option.

•

This option is useful for some variable step methods, as without this limit it is possible for a calculation to run forever, and never get past a certain point in the integration (if the step size were, for example, halving at each step of the numerical integration). This can happen when singularities are present. It is also useful when the function to be evaluated is large and complex or is expensive to evaluate (say, for example, a procedure call that does a numerical sum).



Examples


Consider the DE:
>

DE:=diff(y(t),t)=cos(t^2)*y(t);

${\mathrm{DE}}{\u2254}\frac{{\ⅆ}}{{\ⅆ}{t}}\phantom{\rule[0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{\mathrm{cos}}{}\left({{t}}^{{2}}\right){}{y}{}\left({t}\right)$
 (1) 
A numerical routine for the solution of this DE is obtained by using dsolve and the rkf45 method as follows:
>

ysol:=dsolve({DE, y(0)=1}, numeric, method=rkf45);

${\mathrm{ysol}}{:=}{\mathbf{proc}}\left({\mathrm{x\_rkf45}}\right)\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{\mathbf{end\; proc}}$
 (2) 
The routine is then applied to a value to obtain the numerical solution at that point:
$\left[{t}{=}{5.}{\,}{y}{}\left({t}\right){=}{1.84313059855996}\right]$
 (3) 
Suppose that instead we wanted to integrate out to $t=200$. If we attempt to do this with the generated routine, a problem occurs.
Integrating this DE out to $t=200$ actually requires nearly half a million function evaluations due to the $\mathrm{O}\left({t}^{2}\right)$ increase in the rate of change of the function on the righthand side. This can take some time for computation, so maxfun acts like a stopcheck here.
If we really want this value, we can increase the value of maxfun from its default of $30000$ to $500000$, and obtain a solution as follows:
>

ysol:=dsolve({DE, y(0)=1}, numeric, method=rkf45, maxfun=500000);

${\mathrm{ysol}}{:=}{\mathbf{proc}}\left({\mathrm{x\_rkf45}}\right)\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{\mathbf{end\; proc}}$
 (4) 
$\left[{t}{=}{200.}{\,}{y}{}\left({t}\right){=}{1.87563690274425}\right]$
 (5) 
Alternatively, we could have set maxfun to $0$, and then the routine would proceed until it found the value at the specified point, or recognized a singularity.
The ODE IVP:
>

dsys := {diff(y(t),t)=y(t)^2/(1100*y(t)), y(0)=1/200}:

has a fast blowup singularity just past $t=30$.
If we attempt to numerically integrate this IVP to $t=31$, we get:
>

ysol:=dsolve(dsys, numeric, method=rkf45, abserr=1e7, relerr=1e7);

${\mathrm{ysol}}{:=}{\mathbf{proc}}\left({\mathrm{x\_rkf45}}\right)\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{\mathbf{end\; proc}}$
 (6) 
Or implicitly using the range argument:
>

ysol:=dsolve(dsys, numeric, method=rkf45, abserr=1e7, relerr=1e7,
range=0..31);

${\mathrm{ysol}}{:=}{\mathbf{proc}}\left({\mathrm{x\_rkf45}}\right)\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{\mathbf{end\; proc}}$
 (7) 
This error tells us that integration over this region required more than the default of $30000$ evaluations of the DE. This seems to indicate that we need to increase maxfun. In doing so, we get:
>

ysol:=dsolve(dsys, numeric, method=rkf45, abserr=1e7, relerr=1e7,
maxfun=100000);

${\mathrm{ysol}}{:=}{\mathbf{proc}}\left({\mathrm{x\_rkf45}}\right)\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{\mathbf{end\; proc}}$
 (8) 
we see that we still hit the maxfun barrier, though the integration does proceed a small bit further.
Since there is an actual singularity, this is a case where the singularity detection fails and maxfun set to $30000$ prevents the calculation from running for an excessive time while serving no useful purpose.


