Improved Methods

The default method for numerical solution of IVP problems has been replaced by two new methods written by L. Shampine and R. Corless. The rkf45 method replaces the prior rkf45 code for solution of non-stiff IVP, and the new rosenbrock method is the default method for solution of stiff IVP.

Stopping Conditions

A new option, stop_cond, has been provided for the default methods (rkf45 and rosenbrock). This allows for specification of a method for solving an overdetermined differential algebraic equation problem (DAE), where subsequent failure to satisfy the additional algebraic constraints will halt computation. Examples of constraints or halting conditions include problems where a reactant concentration becomes negative, or where a projectile passes a given boundary in a targeting problem.

The following example is for a projectile targeting problem in rolling hills (sinusoidal):

sys:={diff(x(t),t,t)=0,diff(y(t),t,t)+9.8=0};

$\mathrm{sys}≔\left\{\frac{{\ⅆ}^2}{\ⅆt^2}y\left(t\right)+9.8=0\,\frac{{\ⅆ}^2}{\ⅆt^2}x\left(t\right)=0\right\}$

stopcond:= [y(t)-8*sin(x(t)),x(t)-100];

$\mathrm{stopcond}≔\left[y\left(t\right)-8\sin \left(x\left(t\right)\right)\,x\left(t\right)-100\right]$

IC1 := {x(0)=0,y(0)=10,D(x)(0)=10,D(y)(0)=40};

$\mathrm{IC1}≔\left\{x\left(0\right)=0\,y\left(0\right)=10\,\mathrm{D}\left(x\right)\left(0\right)=10\,\mathrm{D}\left(y\right)\left(0\right)=40\right\}$

ans:=dsolve(sys union IC1, 'numeric', 'stop_cond'=stopcond);

$\mathrm{ans}:=\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

ans(100);

$\left[t=8.25697905978775\,x\left(t\right)=82.5697905978775\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=10.\,y\left(t\right)=6.20841674202050\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=\mathrm{-40.9183947859199}\right]$

IC2 := {x(0)=0,y(0)=10,D(x)(0)=10,D(y)(0)=80};

$\mathrm{IC2}≔\left\{x\left(0\right)=0\,y\left(0\right)=10\,\mathrm{D}\left(x\right)\left(0\right)=10\,\mathrm{D}\left(y\right)\left(0\right)=80\right\}$

ans:=dsolve(sys union IC2, 'numeric', 'stop_cond'=stopcond);

$\mathrm{ans}:=\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

ans(100);

$\left[t=10.0000000000000\,x\left(t\right)=100.000000000000\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=10.\,y\left(t\right)=320.000000000000\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=\mathrm{-18.0000000000000}\right]$

Boundary Value Problems

Maple can numerically solve boundary value problems with the new bvp solver. This solver is designed for non-stiff to mildly stiff BVP and uses grid adaptation and convergence acceleration to obtain the solution. It is also capable of producing high-accuracy solutions using arbitrary precision floating-point arithmetic, and has a method for removable end-point singularities as well as a faster method for problems without singularities.

dsol := dsolve({diff(y(x),x,x)-2*y(x)=0, y(0)=1.2, y(1)=0.9}, 'numeric');

$\mathrm{dsol}:=\mathbf{proc}\left(\mathrm{x\_bvp}\right)...\mathbf{end\; proc}$

dsol(0);

$\left[x=0.\,y\left(x\right)=1.20000000000000\,\frac{\ⅆ}{\ⅆx}y\left(x\right)=\mathrm{-1.25251895272792}\right]$

dsol(0.2);

$\left[x=0.2\,y\left(x\right)=0.994463627112648\,\frac{\ⅆ}{\ⅆx}y\left(x\right)=\mathrm{-0.816528958014159}\right]$

dsol(0.5);

$\left[x=0.5\,y\left(x\right)=0.832942089176703\,\frac{\ⅆ}{\ⅆx}y\left(x\right)=\mathrm{-0.276385195570455}\right]$

dsol(1);

$\left[x=1.\,y\left(x\right)=0.900000000000000\,\frac{\ⅆ}{\ⅆx}y\left(x\right)=0.555701111490340\right]$

High-accuracy solution

Digits:=30:

sys := {diff(y(x),x,x)-2*y(x)=0, y(0)=1.2, y(1)=0.9};

$\mathrm{sys}≔\left\{\frac{{\ⅆ}^2}{\ⅆx^2}y\left(x\right)-2y\left(x\right)=0\,y\left(0\right)=1.2\,y\left(1\right)=0.9\right\}$

exact:=dsolve(sys);

$\mathrm{exact}≔y\left(x\right)=-\frac{3\left(4{\ⅇ}^{-\sqrt{2}}-3\right){\ⅇ}^{\sqrt{2}x}}{10\left(-{\ⅇ}^{-\sqrt{2}}+{\ⅇ}^{\sqrt{2}}\right)}+\frac{3\left(4{\ⅇ}^{\sqrt{2}}-3\right){\ⅇ}^{-\sqrt{2}x}}{10\left(-{\ⅇ}^{-\sqrt{2}}+{\ⅇ}^{\sqrt{2}}\right)}$

v1:=evalf(eval(rhs(exact),x=0.5));

$\mathrm{v1}≔0.832942090833706259044696272732$

dsol := dsolve(sys, 'numeric', 'abserr'=1e-20);

$\mathrm{dsol}:=\mathbf{proc}\left(\mathrm{x\_bvp}\right)...\mathbf{end\; proc}$

v2:=eval(y(x),dsol(0.5));

$\mathrm{v2}≔0.832942090833706259045253674662$

abs(v1-v2);

$5.57401930\times {10}^{\mathrm{-22}}$

Digits:=10:

Speed

The default methods, rkf45 for non-stiff problems and rosenbrock for stiff problems, can be provided an option, $\mathrm{range}=a..b$, to allow precomputation of the ODE solution, allowing interpolation to be used for plotting and obtaining solution values between steps. This results in faster plotting of solutions, higher quality solution plots, better solutions, and avoids prior problems with the dependence of the solution on the resolution of the plot. The solution procedures for these methods store their computations, so if the numerical solution of an IVP problem was calculated for $t=0..100$, any later request for solution values in that range is obtained through interpolation of the already computed solution.

In addition, these methods can run in hardware mode, where the bulk of the computation is being performed by compiled code and evaluation of the ODE system is performed using evalhf. This results in significant speed improvements for problems that are evalhf-compatible.

For example, the following long-time integration of the Lotka-Volterra model is computed in one quarter the time than in prior releases:

dsys := {diff(x(t),t)=x(t)*(1-y(t)),
         diff(y(t),t)=.3*y(t)*(x(t)-1),
         x(0)=1.0, y(0)=0.5};

$\mathrm{dsys}≔\left\{\frac{\ⅆ}{\ⅆt}x\left(t\right)=x\left(t\right)\left(1-y\left(t\right)\right)\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=0.3y\left(t\right)\left(x\left(t\right)-1\right)\,x\left(0\right)=1.0\,y\left(0\right)=0.5\right\}$

dsol := dsolve(dsys, 'numeric', 'maxfun'=0,
               'interpolate'=false);

$\mathrm{dsol}:=\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

dsol(2000);

$\left[t=2000.\,x\left(t\right)=0.356896878813509\,y\left(t\right)=0.657193538022169\right]$

Consistency and Ease of Use

The dsolve/numeric methods have been made more uniform and easier to use.

Nearly all dsolve/numeric methods for IVP understand the $\mathrm{abserr}$ and $\mathrm{relerr}$ options for step error control, the $\mathrm{maxfun}$ option for limiting the amount of work performed by the routine, and the $\mathrm{minstep}$, $\mathrm{maxstep}$, and $\mathrm{initstep}$ options for step size control.

The interface has also been simplified, no longer requiring the names of the dependent variables of the problem (which can be determined from the input), and allowing for specification of the use of a stiff method with the option stiff=true.

deq  := diff(y(t),t,t) + 1001*diff(y(t),t) + 1000*y(t)=0:

ic   := y(0)=1, D(y)(0)=-1:

dsol := dsolve({deq,ic}, 'numeric');

$\mathrm{dsol}:=\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

dsol(50);

Error, (in dsol) cannot evaluate the solution further right of 15.811726, maxfun limit exceeded (see ?dsolve,maxfun for details)

dsol := dsolve({deq,ic}, 'numeric', 'stiff'=true);

$\mathrm{dsol}:=\mathbf{proc}\left(\mathrm{x\_rosenbrock}\right)...\mathbf{end\; proc}$

dsol(50);

$\left[t=50.\,y\left(t\right)=9.37362189971127\times {10}^{\mathrm{-10}}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=\mathrm{-9.37362189971127}\times {10}^{\mathrm{-10}}\right]$

The $\mathrm{mgear}$ method has been removed. The $\mathrm{rosenbrock}$ or $\mathrm{lsode}$ methods should be used instead.

The value=array option is deprecated, as the same functionality is now available as part of the output option.

Precision

The rkf45, rosenbrock, and dverk78 methods can now compute solutions in arbitrary precision Maple floating-point (based on Digits). In prior releases, the taylorseries and gear methods were the only methods capable of arbitrary-precision solutions.

Miscellaneous

The current IVP methods (with the exception of taylorseries) can now handle complex-valued problems.

dsol := dsolve({diff(y(t),t) + I*y(t) = 0, y(0)=I},
               'numeric'):

dsol(Pi/4);

$\left[t=0.785398163397448\,y\left(t\right)=0.707106875083858+0.707106902685889\mathrm{I}\right]$

dsol := dsolve({diff(y(t),t) + I*y(t) = 0, y(0)=I},
               'numeric', 'method'='rosenbrock'):

dsol(Pi/4);

$\left[t=0.785398163397448\,y\left(t\right)=0.707106735017522+0.707106826931253\mathrm{I}\right]$

dsol := dsolve({diff(y(t),t) + I*y(t) = 0, y(0)=I},
               'numeric', 'method'='gear'):

dsol(Pi/4);

$\left[t=0.785398163397448\,y\left(t\right)=0.707106781179390+0.707106781190284\mathrm{I}\right]$

dsol := dsolve({diff(y(t),t) + I*y(t) = 0, y(0)=I},
               'numeric', 'method'='dverk78'):

dsol(Pi/4);

$\left[t=0.785398163397448\,y\left(t\right)=0.707106781184395+0.707106781187889\mathrm{I}\right]$

dsol := dsolve({diff(y(t),t) + I*y(t) = 0, y(0)=I},
               'numeric', 'method'='lsode'):

dsol(Pi/4);

$\left[t=0.785398163397448\,y\left(t\right)=0.707107828482804+0.707106783632396\mathrm{I}\right]$

The solution outputs when using the $\mathrm{output}=\mathrm{listprocedure}$ option are now more flexible, and can be used directly with other Maple commands, such as fsolve and plot.

dsol := dsolve({diff(y(t),t,t)+y(t)=0,y(0)=1,D(y)(0)=0},
               'numeric', 'abserr'=1e-9, 'relerr'=1e-9,
               'output'=listprocedure);

$\mathrm{dsol}≔\left[t=\mathbf{proc}\left(t\right)...\mathbf{end\; proc}\,y\left(t\right)=\mathbf{proc}\left(t\right)...\mathbf{end\; proc}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=\mathbf{proc}\left(t\right)...\mathbf{end\; proc}\right]$

ysol := eval(y(t),dsol);

$\mathrm{ysol}:=\mathbf{proc}\left(t\right)...\mathbf{end\; proc}$

fsolve(ysol(t)=0,t,1..2);

$1.570796327$

evalf(Pi/2);

$1.570796327$

A long-standing bug, relating to integration of IVP on both sides of the initial condition has been fixed for all methods.

The most common forms for invoking numerical integration remain unchanged, namely

evalf(Int(expr, x=a..b))

`evalf/int`(expr, x=a..b))

evalf(int(expr, x=a..b))

where the last form only invokes numerical integration if the invocation of symbolic integration returns an unevaluated integral. There is a variation of the above forms in the case where the first argument (the integrand) is specified as a procedure or a Maple operator, in which case the second argument must be the range a..b and not an equation.

The syntax for specifying additional options has been modified in Maple 7 to accept equation forms. (For backward compatibility, some options are accepted as values rather than equations, as specified below.) One or more of the following options may be specified as additional arguments:

The specification method = <name> (or simply <name>) indicates a particular numerical integration method to be applied. For the list of methods which may be specified see evalf,int. By default, a hybrid symbolic-numeric strategy is applied.

The specification digits = <posint> (or simply <posint>) indicates the number of digits of precision for the computation. Some additional guard digits are carried during the computation to attempt to achieve a result with <posint> correct digits (although a larger tolerance can be specified by using the 'epsilon' option). By default, the Maple environment variable Digits specifies the precision for the computation.

The specification epsilon = <numeric> specifies the relative error tolerance for the computed result. The routines attempt to achieve a final result with a relative error less than this value. By default, the relative error tolerance which the routines attempt to achieve for the final result is

where digits is the precision specified for the computation. In attempting to achieve this accuracy, the working value of Digits is increased as deemed necessary. It is an error to specify 'epsilon' smaller than the default value above.

In Maple 7, compiled routines from the NAG (Numerical Algorithms Group) subroutine library are automatically invoked to perform quadrature (numerical integration) tasks, when appropriate. This results in significant speed improvements compared with earlier versions of Maple. (See below for some timing information.)

In the default case (no particular method specified), the Maple integration problem is first passed to NAG integration routines if Digits is not too large (i.e., if $\mathrm{Digits}\le \mathrm{evalhf}\left(\mathrm{Digits}\right)$). The NAG routines are in a compiled C library and hence operate at hardware floating-point speed. If the NAG routines cannot perform the integration, then some singularity handling may be performed and control may pass back to the NAG routines with a modified problem. Native Maple routines are invoked if the NAG routines cannot solve the problem (e.g., for larger values of Digits and for integrands involving functions unknown to NAG).

Various modifications in the hybrid symbolic-numeric solution strategy for integrands with singularities have been implemented, driven by a study of cases which were not handled (or were inefficiently handled) in previous versions. Thus a wider class of integrands can now be integrated successfully.

For example, integrands containing various nonanalytic functions are treated by breaking the interval of integration at the points at which the function is not analytic and by separating real and imaginary parts of the integrand. The following nonanalytic functions known to Maple are treated: piecewise, Dirac, Heaviside, abs, signum, csgn, trunc, frac, floor, ceil, round, ilog10, argument .

The following two examples will illustrate. By setting ${\mathrm{infolevel}}_{\mathrm{evalf/int}}≔2$ one can see how the interval of integration is split at each potential singularity of the integrand $\mathrm{f1}$ and the integrand $\mathrm{f2}$ is split into its real and imaginary parts.

infolevel[`evalf/int`] := 2:

f1 := ceil(cos(x^2-1));

$\mathrm{f1}≔⌈\cos \left(x^2-1\right)⌉$

evalf(Int(f1, x=0..Pi));

evalf/int/control:   integrating on 0 .. Pi the integrand

evalf/int/control:   non-analytic function detectedevalf/int/pieces:   attempt to break up the interval of integration
evalf/int/control:   integrating on 0 .. 1 the integrand

evalf/int/control:   integrating on 1 .. 1/2*(2*Pi+4)^(1/2) the integrand

evalf/int/control:   integrating on 1/2*(2*Pi+4)^(1/2) .. (1+Pi)^(1/2) the integrand

evalf/int/control:   integrating on (1+Pi)^(1/2) .. 1/2*(6*Pi+4)^(1/2) the integrand

evalf/int/control:   integrating on 1/2*(6*Pi+4)^(1/2) .. (2*Pi+1)^(1/2) the integrand

evalf/int/control:   integrating on (2*Pi+1)^(1/2) .. 1/2*(10*Pi+4)^(1/2) the integrand

evalf/int/control:   integrating on 1/2*(10*Pi+4)^(1/2) .. Pi the integrand

$⌈\cos \left(x^2-1\right)⌉$

$1$

$1$

$0$

$0$

$1$

$1$

$0$

$2.188873934$

f2 := conjugate(sin(x+I*arctan(x^2,sin(x))));

$\mathrm{f2}≔\sin \left(\stackrel{\&conjugate0;}{\left(x+\mathrm{I}\arctan \left(x^2\&comma;\sin \left(x\right)\right)\right)}\right)$

evalf(Int(f2, x=0..2*Pi));

evalf/int/control:   integrating on 0 .. 2*Pi the integrand

evalf/int/control:   attempt to split integrand into Re and Im partsevalf/int/control:   integrating on 0 .. 2*Pi the integrand

Control:   Entering NAGIntControl:   trying d01ajc (nag_1d_quad_gen)
d01ajc:   "trying evalhf callbacks"
d01ajc:   result=-1.65939040794013204
d01ajc:   abserr=.204184750295993220e-9; num_subint=2; fun_count=63
Control:   result=-1.65939040794013204
evalf/int/control:   integrating on 0 .. 2*Pi the integrand

Control:   Entering NAGIntControl:   trying d01ajc (nag_1d_quad_gen)
d01ajc:   "trying evalhf callbacks"
d01ajc:   result=1.55748944065209516
d01ajc:   abserr=.728975743687537407e-10; num_subint=2; fun_count=63
Control:   result=1.55748944065209516

$\sin \left(\stackrel{\&conjugate0;}{\left(x+\mathrm{I}\arctan \left(x^2\&comma;\sin \left(x\right)\right)\right)}\right)$

$\sin \left(x\right)\cosh \left(\arctan \left(x^2\&comma;\sin \left(x\right)\right)\right)$

$-\cos \left(x\right)\sinh \left(\arctan \left(x^2\&comma;\sin \left(x\right)\right)\right)$

$\mathrm{-1.659390408}+1.557489441\mathrm{I}$

The following two timing tests are examples of the timing improvements for numerical integration which have been achieved in Maple 7 due to the incorporation of compiled NAG routines. For integration problems on a finite interval with no integrand singularities, the speedup is typically a factor of 20 to 30.

ex1

evalf(Int(ln(2*x^2)*x^4 / (1 - exp(2*x^2)), x = -infinity..infinity));

evalf/int/control:   integrating on -infinity .. infinity the integrand

evalf/int/control:   "attempting Levin method"Control:   Entering NAGInt
trying d01amc (nag_1d_quad_inf)
d01amc:   "trying evalhf callbacks"
d01amc:   result=-.653212119735574331e-1
d01amc:   abserr=.222032947583272744e-10; num_subint=12; fun_count=345
result=-.653212119735574331e-1
Control:   Entering NAGInt
trying d01amc (nag_1d_quad_inf)
d01amc:   "trying evalhf callbacks"
d01amc:   result=-.653212119735574331e-1
d01amc:   abserr=.222032947583272744e-10; num_subint=12; fun_count=345
result=-.653212119735574331e-1

$\frac{\ln \left(2x^2\right)x^4}{1-{\&ExponentialE;}^{2x^2}}$

$\mathrm{-0.1306424239}$

ex2

evalf(Int(1/exp(x*ln(x)), x=1..infinity));

evalf/int/control:   integrating on 1 .. infinity the integrand

evalf/int/control:   "attempting Levin method"evalf/int/IsolateInt:   splitting interval based on points: [1, 10.8462559396596, infinity]
evalf/int/control:   integrating on 10.8462559396596 .. infinity the integrand

evalf/int/control:   "attempting Levin method"Control:   Entering NAGInt
trying d01amc (nag_1d_quad_inf)
d01amc:   "trying evalhf callbacks"
d01amc:   result=.173116322605122934e-11
d01amc:   abserr=.860731693119929988e-14; num_subint=1; fun_count=15
result=.173116322605122934e-11
evalf/int/control:   integrating on 1 .. 10.8462559396596 the integrand

Control:   Entering NAGIntControl:   trying d01ajc (nag_1d_quad_gen)
d01ajc:   "trying evalhf callbacks"
d01ajc:   result=.704169960435743203
d01ajc:   abserr=.472390388719768988e-10; num_subint=2; fun_count=63
Control:   result=.704169960435743203

$\frac{1}{{\&ExponentialE;}^{x\ln \left(x\right)}}$

${\left(\frac{1}{x}\right)}^x$

$\frac{1}{{\&ExponentialE;}^{x\ln \left(x\right)}}$

$0.7041699604$