Error Control for Numerical Solution of IVPs in Maple
Description
Examples
dsolve[numeric] approximates the solution of an initial value problem (IVP) for a system of ordinary differential equations (ODEs). Although ODEs arise in various forms, they can always be rewritten as a first order system, so to explain error control, it is assumed that the ODEs have the form
ⅆⅆxyx=fx,yx,ya=y0
for x in the range a..b. Here x is the independent variable and yx is the vector of dependent variables.
Numerical methods solve this IVP by approximating the solution successively at a discrete set of points x0=a<x1<x2<...<xm=b. Using an approximation yi to yxi, and possibly approximations at a few xj prior to xi, an approximation yi+1 is computed to yxi+1. This is called taking a step to xi+1 with a step of size h=xi+1−xi.
At each step the solver makes a truncation and/or discretization error that depends on the method and the length of the step. The cumulative effect of these errors depends on the stability of the IVP near the solution yx. If the IVP is stable near the solution (solutions with nearby initial data do not diverge from one another), errors are not amplified, but if the IVP is unstable near the solution (solutions with nearby initial data diverge), errors are amplified.
Adaptive solvers estimate the discretization error at each step and control it by adjusting the step size. The cumulative error in the numerical solution depends on the stability of the IVP, but for moderately stable problems, the solvers are tuned so that the error is comparable to the tolerances on the discretization error. For the efficient solution of an IVP, the solvers use the largest step size that results in a discretization error smaller than the tolerances. That is, a large step size is used if the solution is easy to approximate and a small one if it is difficult.
Unfortunately, there is no "best" numerical method for IVPs. A method for which the discretization error is Ohp is said to be of order p. If the step size h is sufficiently small, a method of higher order has a smaller discretization error. Because methods of higher order are more expensive, this is an advantage only if the tolerances are sufficiently small. For greater accuracy, use a higher order method. In particular, dsolve[rkf45] is generally more efficient than the higher order dsolve[dverk78] for modest tolerances (the default), but dsolve[dverk78] is generally more efficient for stringent tolerances. The dsolve[lsode] code adapts the order, as well as the step size, so as to be efficient over a wide range of tolerances.
The Maple numerical IVP solvers control the discretization error by means of the options abserr, relerr, minstep, maxstep, and initstep. Note that the classical methods do not estimate the discretization error or vary the step size to control it.
The options that pertain to error control are:
abserr
relerr
initstep
maxstep
minstep
Not all options apply to all methods, and exceptions are noted in the discussion of each option. In all cases, the default values are specific to each method. For additional information, consult the individual method help pages.
abserr and relerr are tolerances for the discretization error. abserr is an absolute error tolerance and relerr is a relative error tolerance. The exact meaning of these tolerances is dependent upon the solver. To explain this, let the dependent variables at the i-th step be yij,j=1..n and the estimated discretization error be approxerrj.
For the rkf45, ck45, rosenbrock, the related DAE extension methods, the taylorseries method, and the mebdfi method, the inequality
approxerrj≤abserr+relerryij
must be satisfied for all yij simultaneously, j=1..n.
Note: some of the solvers (all but taylorseries mentioned above) allow for a per-component absolute error tolerance, so in this case the inequality becomes:
approxerrj≤abserrj+relerryij
For the gear and dverk78 methods, the inequality
approxerrj≤maxabserr,relerryij
lsode measures the error in the sense of root-mean-square:
∑j=1napproxerrj2abserr+relerryij2n≤1
The other options, minstep, maxstep, and initstep, are generally intended as fine tuning options and must be used with care.
The minstep parameter places a minimum on the size of the steps taken by the solver, and forces an error if the solver cannot achieve the required error tolerance without reducing the step below the minimum. This can be used to recognize singularities, and can also be used to limit the cost of the of the computation, though a better way to accomplish this is to limit the number of evaluations of fx,y, see dsolve[maxfun].
The maxstep parameter places a maximum on the size of the steps taken by the solver, so that even if the step control indicates that a larger step is possible, the size of the step will not exceed the imposed limit. This can be used to ensure that the solver does not lose the scale of the problem.
You can specify to the solver the scale of the problem near the initial point x=a by supplying an initial step size as initstep. The solver uses this step size if the error of the step is acceptable. Otherwise, it reduces the step size and tries again.
The minstep and maxstep options are not available for the solvers rkf45, ck45, rosenbrock, and the related DAE extension methods, as they require control of these parameters to provide a good solution.
dsys≔diffyx,x,x=−yx,y0=0,Dy0=1
dsys≔ⅆ2ⅆx2yx=−yx,y0=0,Dy0=1
The first example solves this IVP with the default rkf45 method and default error tolerances. The errf function measures the difference between the computed solution and the exact solution. This example shows how the error can vary over the interval of integration. The remaining examples show the effect of tighter tolerances on the error.
The rkf45 method with default error tolerances. Note that we need the option maxfun=0 to remove the limit on the number of function evaluations performed for the integration out to t=900:
dsol1≔dsolvedsys,numeric,maxfun=0
dsol1 ≔ procx_rkf45...end proc
errf≔x↦evalfabssinx−rhsdsol1x2
errf≔x↦evalfsinx−rhsdsol1x2
errf9,errf90,errf900
4.46240948692722×10−7,8.91718218909432×10−6,0.0000929084983034567
The following plot shows the growth of the global error over the first portion of the interval of computation.
ploterrf,0..9
The rkf45 method with tighter error tolerances:
dsol2≔dsolvedsys,numeric,abserr=1.×10−8,relerr=1.×10−8,maxfun=0
dsol2 ≔ procx_rkf45...end proc
errf≔x↦evalfabssinx−rhsdsol2x2
errf≔x↦evalfsinx−rhsdsol2x2
1.15294584435155×10−8,2.42536584593722×10−7,2.62007740259307×10−6
The dverk78 method with the same tolerances. This solver is more efficient than rkf45 at these tolerances, so it is not necessary to increase the default maxfun.
dsol3≔subsdsolvedsys,numeric,output=listprocedure,method=dverk78,abserr=1.×10−8,relerr=1.×10−8,yx
dsol3 ≔ procx...end proc
errf≔x↦evalfabssinx−dsol3x
errf≔x↦evalfsinx−dsol3x
4.93375695853615×10−9,1.66151348235388×10−9,4.79303242650886×10−7
This solution is sufficiently accurate that, to compute the error accurately, we must use more than the default 10 digits of accuracy in Maple.
Digits≔truncevalhfDigits
The dverk78 method with even tighter error tolerances:
dsol4≔subsdsolvedsys,numeric,output=listprocedure,method=dverk78,abserr=1.×10−10,relerr=1.×10−10,yx
dsol4 ≔ procx...end proc
errf≔x↦evalfabssinx−dsol4x
errf≔x↦evalfsinx−dsol4x
1.33253463818761×10−10,2.34880781491142×10−10,4.41993575073241×10−9
See Also
dsolve[ck45]
dsolve[classical]
dsolve[dverk78]
dsolve[gear]
dsolve[lsode]
dsolve[maxfun]
dsolve[numeric,IVP]
dsolve[numeric]
dsolve[rkf45]
dsolve[rosenbrock]
dsolve[taylorseries]
Download Help Document
