|
Option maxfun for Numerical Solution of ODE Initial Value Problems
|
|
Description
|
|
|
•
|
Invoking the dsolve function with the numeric option and  (default),  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:
|
|
•
|
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:
|
>
|
|
|
| (1) |
A numerical routine for the solution of this DE is obtained by using dsolve and the rkf45 method as follows:
|
>
|
|
|
| (2) |
The routine is then applied to a value to obtain the numerical solution at that point:
|
>
|
|
|
![[t = 5., y(t) = 1.84313059855996]](/support/helpjp/helpview.aspx?si=3889/file06712/math164.png)
| (3) |
Suppose that instead we wanted to integrate out to  . If we attempt to do this with the generated routine, a problem occurs.
|
>
|
|
Integrating this DE out to  actually requires nearly half a million function evaluations due to the  increase in the rate of change of the function on the right-hand 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  to  , and obtain a solution as follows:
|
>
|
|
|
| (4) |
|
>
|
|
|
![[t = 200., y(t) = 1.87563259797483]](/support/helpjp/helpview.aspx?si=3889/file06712/math208.png)
| (5) |
Alternatively, we could have set maxfun to  , and then the routine would proceed until it found the value at the specified point, or recognized a singularity.
The ODE IVP:
|
>
|
|
has a fast blow-up singularity just past  .
If we attempt to numerically integrate this IVP to  , we get:
|
>
|
|
|
| (6) |
|
>
|
|
Or implicitly using the range argument:
|
>
|
|
|
Warning, cannot evaluate the solution further right of 30.685302, maxfun limit exceeded (see ?dsolve,maxfun for details)
| |
|
| (7) |
This error tells us that integration over this region required more than the default of  evaluations of the DE. This seems to indicate that we need to increase maxfun. In doing so, we get:
|
>
|
|
|
| (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  prevents the calculation from running for an excessive time while serving no useful purpose.
|
|
