Documentation |
Numerical solution of an ordinary differential equation
This functionality does not run in MATLAB.
numeric::odesolve(f, t_{0} .. t, Y_{0}, <method>, <RelativeError = rtol>, <AbsoluteError = atol>, <Stepsize = h>, <MaxStepsize = h_{max}>, <Alldata = n>, <Symbolic>) numeric::odesolve(t_{0} .. t, f, Y_{0}, <method>, <RelativeError = rtol>, <AbsoluteError = atol>, <Stepsize = h>, <MaxStepsize = h_{max}>, <Alldata = n>, <Symbolic>)
numeric::odesolve(f, t_{0}..t, Y_{0}) returns a numerical approximation of the solution Y(t) of the first order differential equation (dynamical system) , Y(t_{0}) = Y_{0} with and .
numeric::odesolve is a general purpose solver able to deal with initial value problems of various kinds of ordinary differential equations. Equations of order p can be solved by numeric::odesolve after reformulation to dynamical system form. This can always be achieved by writing the equation as a first order system
for the vector . See Example 4.
The following single-step Runge-Kutta-type methods are implemented:
EULER1 (order 1)
RKF43 (order 3)
xRKF43 (order 3)
RKF34 (order 4)
xRKF34 (order 4)
RK4 (order 4)
RKF54a (order 4)
RKF54b (order 4)
DOPRI54 (order 4)
xDOPRI54 (order 4)
CK54 (order 4)
xRKF54a (order 4)
xRKF54b (order 4)
xCK54 (order 4)
RKF45a (order 5)
RKF45b (order 5)
DOPRI45 (order 5)
CK45 (order 5)
xRKF45a (order 5)
xRKF45b (order 5)
xDOPRI45 (order 5)
xCK45 (order 5)
DOPRI65 (order 5)
xDOPRI65 (order 5)
DOPRI56 (order 6)
xDOPRI56 (order 6)
BUTCHER6 (order 6)
RKF87 (order 7)
xRKF87 (order 7)
DOPRI87 (order 7)
xDOPRI87 (order 7)
RKF78 (order 8)
xRKF78 (order 8)
DOPRI78 (order 8)
xDOPRI78 (order 8)
GAUSS(s) (order 2 s)
GAUSS = s
For the Gauss methods, GAUSS(s) is equivalent to GAUSS = s. The positive integer s indicates the number of stages. The order of the s stage Gauss method is 2 s.
The utility function numeric::ode2vectorfield may be used to produce the input parameters f, t_{0}, Y_{0} from a set of differential expressions representing the ODE. See Example 1.
The input data t_{0}, t and Y_{0} must not contain symbolic objects which cannot be converted to floating point values via float. Numerical expressions such as , etc. are accepted.
The vector field f defining the dynamical system must be represented by a procedure with two input parameters: the scalar time t and the vector Y. numeric::odesolve internally calls this function with real floating-point values t and a list Y of floating-point values. It has to return the vector f(t, Y) either as a list or as a 1-dimensional array. The output of f may contain numerical expressions such as π, etc. However, all values must be convertible to real or complex floating point numbers by float.
Autonomous systems, where f(t, Y) does not depend on t, must also be represented by a procedure with 2 arguments t and Y.
Scalar functions Y also must be represented by a list or an array with one element. For instance, the input data for the scalar initial value problem may be of the form
f := proc(t,Y) | /* Y is a 1-dimensional vector */ |
local y; | /* represented by a list with */ |
begin | /* one element: Y = [y]. */ |
y := Y[1]; | |
[t*sin(y)] | /* the output is a list with 1 element */ |
end_proc: | |
Y0 := [1]: | /* the initial value */ |
The numerical precision is controlled by the global variable DIGITS: an adaptive control of the step size keeps local relative discretization errors below rtol=10^-DIGITS, unless a different tolerance is specified via the option RelativeError = rtol. For small values of the solution vector Y, the absolute discretization error can be bounded by the threshold atol specified via the option AbsoluteError = atol.
If AbsoluteError is not specified, only relative discretization errors are controlled and kept below rtol.
The error control may be switched off by specifying a fixed Stepsize = h.
Note: Only local errors are controlled by the adaptive mechanism. No control of the global error is provided! |
With Y := t -> numeric::odesolve(f, t_0..t, Y_0), the numerical solution can be repesented by a MuPAD^{®} function: the call Y(t) will start the numerical integration from t_{0} to t. A more sophisticated form of this function may be generated via Y := numeric::odesolve2(f, t_{0}, Y_{0}).
This equips Y with a remember mechanism that uses previously computed values to speed up the computation. See Example 2.
For systems of the special form with a matrix valued function f(t, Y), there is a special solver numeric::odesolveGeometric which preserves geometric features of the system more faithfully than numeric::odesolve.
The function is sensitive to the environment variable DIGITS, which determines the numerical working precision.
We compute the numerical solution y(10) of the initial value problem , y(0) = 2:
f := proc(t, Y) begin [t*sin(Y[1])] end_proc: numeric::odesolve(f, 0..10, [2])
Alternatively, the utility function numeric::ode2vectorfield can be used to generate the input parameters in a more intuitive way:
[f, t0, Y0] := [numeric::ode2vectorfield({y'(t) = t*sin(y(t)), y(0) = 2}, [y(t)])]
numeric::odesolve(f, t0..10, Y0)
delete f, t0, Y0:
We consider , y(0) = 1. The numerical solution may be represented by the function
Y := t -> numeric::odesolve((t,Y) -> Y, 0..t, [1]):
Calling Y(t) starts the numerical integration:
Y(-5), Y(0), Y(1), Y(PI)
delete Y:
We compute the numerical solution Y(π) = [x(π), y(π)] of the system
.
f := (t, Y) -> [Y[1] + Y[2], Y[1] - Y[2]]: Y0 := [1, I]: numeric::odesolve(f, 0..PI, Y0)
The solution of a linear dynamical system with a constant matrix A is . The solution of the system above can also be computed by:
t := PI: tA := array(1..2, 1..2, [[t, t], [t, -t]]): numeric::expMatrix(tA, Y0)
delete f, Y0, t, tA:
We compute the numerical solution y(1) of the differential equation with initial conditions y(0) = 0, . The second order equation is converted to a first order system for the vector :
.
f := proc(t, Y) begin [Y[2], Y[1]^2] end_proc: Y0 := [0, 1]: numeric::odesolve(f, 0..1, Y0)
delete f, Y0:
We demonstrate how numerical data can be obtained on a user defined time mesh t[i]. The initial value problem is , y(0) = 1, the sample points are t_{0} = 0.0, t_{1} = 0.1, …, t_{100} = 10.0. First, we define the differential equation and the initial condition:
f := (t, Y) -> [sin(t) - Y[1]]: Y[0] := [1]:
We define the time mesh:
for i from 0 to 100 do t[i] := i/10 end_for:
The equation is integrated iteratively from t[i-1] to t[i] with a working precision of 4 significant decimal places:
DIGITS := 4: for i from 1 to 100 do Y[i] := numeric::odesolve(f, t[i-1]..t[i], Y[i-1]) end_for:
The following mesh data are produced:
[t[i], Y[i]] $ i = 0..100
These data can be displayed by a list plot:
plotpoints := [[t[i], op(Y[i])] $ i = 0..100]: plot(plot::PointList2d(plotpoints, PointColor = RGB::Black)):
The same plot can be obtained directly via plot::Ode2d:
plot(plot::Ode2d( [t[i] $ i = 0..100], f, Y[0], [(t, Y) -> [t, Y[1]], Style = Points, Color = RGB::Black]))
delete f, t, DIGITS, Y, plotpoints:
We compute the numerical solution y(1) of by the classical 4-th order Runge-Kutta method RK4. By internal local extrapolation, its effective order is 5:
f := (t, Y) -> Y: DIGITS := 13: numeric::odesolve(f, 0..1, [1], RK4)
Next, we use local extrapolation xRKF78 of the 8-th order submethod of the Runge-Kutta-Fehlberg pair RKF78. This scheme has effective order 9:
numeric::odesolve(f, 0..1, [1], xRKF78)
Both methods yield the same result because of the internal adaptive error control. However, due to its higher order, the method xRKF78 is faster.
delete f, DIGITS:
We consider the stiff ODE . The default method DOPRI78 is explicit and not very efficient in solving very stiff problems:
f := (t, Y) -> [10^4*(cos(t) - Y[1])]: t0 := time(): numeric::odesolve(f, 0..1, [1]), (time() - t0)*msec
We use the implicit A-stable method GAUSS(6). For this stiff problem, it is more efficient than the default method DOPRI78:
t0 := time(): numeric::odesolve(f, 0..1, [1], GAUSS(6)), (time() - t0)*msec
delete t0:
We consider the initial value problem , y(0) = 1. We note that the numerical evaluation of the right hand side of the equation suffers from cancellation effects, when |y| is small.
f := (t, Y) -> [-10^20*Y[1]*(1 - cos(Y[1]))]: Y0 := [1]:
We first attempt to compute y(1) with a working precision of 6 digits using the default setting RelativeError = 10^-DIGITS=10^(-6):
DIGITS := 6: numeric::odesolve(f, 0..1, Y0)
Due to numerical round-off in the internal steps, no digit of this result is correct. Next, we use a working precision of 20 significant decimal places to eliminate roundoff effects:
DIGITS := 20: numeric::odesolve(f, 0..1, Y0, RelativeError = 10^(-6))
Since relative local discretization errors are of the magnitude 10^{- 6}, not all displayed digits are trustworthy. We finally use a working precision of 20 digits and constrain the local relative discretization errors by the tolerance 10^{- 10}:
numeric::odesolve(f, 0..1, Y0, RelativeError = 10^(-10))
delete f, Y0, DIGITS:
We compute the numerical solution y(1) of , y(0) = 1 with various methods and various constant step sizes. We compare the result with the exact solution .
f := (t, Y) -> Y: Y0 := [1]:
We first use the Euler method of order 1 with two different step sizes:
Y := numeric::odesolve(f, 0..1, Y0, EULER1, Stepsize = 0.1): Y, globalerror = float(exp(1)) - Y[1]
Decreasing the step size by a factor of 10 should reduce the global error by a factor of 10. Indeed:
Y := numeric::odesolve(f, 0..1, Y0, EULER1, Stepsize = 0.01): Y, globalerror = float(exp(1)) - Y[1]
Next, we use the classical Runge-Kutta method of order 4 with two different step sizes:
Y := numeric::odesolve(f, 0..1, Y0, RK4, Stepsize = 0.1): Y, globalerror = float(exp(1)) - Y[1]
Decreasing the step size by a factor of 10 in a 4-th order scheme should reduce the global error by a factor of 10^{4}. Indeed:
Y := numeric::odesolve(f, 0..1, Y0, RK4, Stepsize = 0.01): Y, globalerror = float(exp(1)) - Y[1]
delete f, Y0, Y:
We integrate , y(0) = 1 over the interval t ∈ [0, 0.99] with a working precision of 4 digits. The exact solution is . Note the singularity at t = 1.
DIGITS := 4: f := (t, Y) -> [Y[1]^2]: Y0 := [1]:
The option Alldata, equivalent to Alldata = 1, yields all mesh data generated during the internal adaptive process:
numeric::odesolve(f, 0..0.99, Y0, Alldata)
With Alldata = 2, only each second point is returned:
numeric::odesolve(f, 0..0.99, Y0, Alldata = 2)
One can control the time mesh using the option Stepsize = h:
numeric::odesolve(f, 0..0.99, Y0, Stepsize=0.1, Alldata = 1)
However, with the option Stepsize = h, no automatic error control is provided by numeric::odesolve. Note the poor approximation y(t) = 94.3 for t = 0.99 (the exact value is y(0.99) = 100.0). The next computation with smaller step size yields better results:
numeric::odesolve(f, 0..0.99, Y0, Stepsize = 0.01, Alldata = 10)
Example 5 demonstrates how accurate numerical data on a user defined time mesh can be generated using the automatic error control of numeric::odesolve.
delete DIGITS, f, Y0:
The second order equation is written as the dynamical system , for the vector Y = [y, z]. A single symbolic step
of the Euler method is computed:
f := proc(t, Y) begin [Y[2], -sin(Y[1])] end_proc: numeric::odesolve(f, t0..t0+h, [y0, z0], EULER1, Symbolic)
delete f:
BUTCHER6, CK45, CK54, DOPRI45, DOPRI54, DOPRI56, DOPRI65, DOPRI78, DOPRI87, EULER1, RK4, RKF34, RKF43, RKF45a, RKF45b, RKF54a, RKF54b, RKF78, RKF87, xCK45, xCK54, xDOPRI45, xDOPRI54, xDOPRI56, xDOPRI65, xDOPRI78, xDOPRI87, xRKF34, xRKF43, xRKF45a, xRKF45b, xRKF54a, xRKF54b, xRKF78, xRKF87 |
Name of the Runge-Kutta scheme. See Example 6. For details on these schemes, see the Algorithms section. | |
GAUSS |
Name of the Runge-Kutta scheme specified as GAUSS(s) or GAUSS = s. The methods GAUSS(s) or, equivalently, GAUSS = s are the implicit Gauss methods with s stages of order 2 s. These methods are implicit A-stable schemes. The time steps are rather costly to compute. The Gauss methods are useful for integrating stiff ODEs. For non-stiff ODEs, there is usually no need to change the default method DOPRI78. This method is an embedded Runge-Kutta pair of orders 7 and 8. Further, the Gauss methods are symplectic methods. When used with constant step size (Stepsize = h), numerical integration of Hamiltonian systems benefits from this property. See Example 7. | |
RelativeError, AbsoluteError |
Option specified as RelativeError = rtol forces internal numerical Runge-Kutta steps to use step sizes with relative local discretization errors below rtol. This tolerance must be a positive numerical real value not smaller than . The default tolerance is RelativeError = 10^(-DIGITS). Option specified as AbsoluteError = atol forces internal numerical Runge-Kutta steps to use step sizes with absolute local discretization errors below atol. This tolerance must be a nonnegative numerical real value. The default tolerance is AbsoluteError = 10^(-10*DIGITS). The internal control mechanism estimates the local discretization error of a Runge-Kutta step and adjusts the step size adaptively to keep this error below the specified tolerances rtol or atol, respectively. The code uses the criterion
For accepting a solution vector Y. Roughly speaking, the relative error is controlled when the solution Y is sufficiently large. For very small solution values Y, absolute discretization errors are kept below the threshold atol. Specify AbsoluteError = 0 if only control of the relative discretization errors is desired. The error control may be switched off by specifying a fixed Stepsize = h. The default setting of , ensures that the local discretization errors are of the same order of magnitude as numerical roundoff. Usually there is no need to use these options to change this setting. However, occasionally the numerical evaluation of the Runge-Kutta steps may be ill-conditioned or step sizes are so small that the time parameter cannot be incremented by the step size within working precision. In such a case these options may be used to bound the local discretization errors and use a higher working precision given by DIGITS. Only positive real numerical values are accepted.
See Example 8. | |
Stepsize |
Option, specified as Stepsize = h Switches off the internal error control and uses a Runge-Kutta iteration with constant step size h which must be a positive numerical value. By default, numeric::odesolve uses an adaptive step size control mechanism in the numerical iteration. The option Stepsize = h switches off this adaptive mechanism and uses the Runge-Kutta method specified (or the default method DOPRI78) with constant step size h. A final step with smaller step size is used to match the end t of the integration interval t_0..t if is not an integer.
There is usually no need to invoke this option. However, occasionally the built-in adaptive error control mechanism may fail when integrating close to a singularity. In such a case this option may be used to customize a control mechanism for the global error by using different step sizes and investigating the convergence of the corresponding results. Cf. Example 9. | |
MaxStepsize |
Option, specified as MaxStepsize = hmax Restricts adaptive step sizes to values not larger than h_max; h_max must be a positive numerical value. By default, numeric::odesolve uses an adaptive step size control mechanism in the numerical iteration. The option MaxStepsize = h_{max} restricts the adaptive step size to values no larger than h_{max}. If a larger stepsize h is requested explicitly by Stepsize = h, the option MaxStepsize = h_{max} reduces h to h_{max}. | |
Alldata |
Option, specified as Alldata = n Makes numeric::odesolve return a list of numerical mesh points generated by the internal Runge-Kutta iteration. The integer n controls the size of the output list. With this option, numeric::odesolve returns a list of numerical mesh points [[t_{0}, Y_{0}], [t_{1}, Y_{1}], …, [t, Y(t)]] generated by the internal Runge-Kutta iteration. The integer n controls the size of the output list. For n = 1, all internal mesh points are returned. This case may also be invoked by entering the simplified option Alldata, which is equivalent to Alldata = 1. For n > 1, only each n-th mesh point is stored in the list. The list always contains the initial point [t_{0}, Y_{0}] and the final point [t, Y(t)]. For n ≤ 0, only the data [[t_{0}, Y_{0}], [t, Y(t)]] are returned. The output list may be useful to inspect the internal numerical process. Also further graphical processing of the mesh data may be useful. Cf. Example 10. | |
Symbolic |
Makes numeric::odesolve return a vector of symbolic expressions representing a single symbolic step of the Runge-Kutta iteration. The call numeric::odesolve(f, t_{0}..t, Y_{0}, < method >, Symbolic) returns a vector (list or array) of expressions representing a single step of the numerical scheme with step size t - t_{0}. In this mode symbolic values for t_{0}, t and the components of Y_{0} are accepted. This option may be useful if the specified numerical method applied to a given differential equation is to be investigated symbolically. Cf. Example 11. |
The solution vector Y(t) is returned as a list or as a 1-dimensional array of floating-point values. The type of the result vector coincides with the type of the input vector Y_{0}.
With the option Alldata, a list of mesh data is returned.
All methods presently implemented are adaptive versions of Runge-Kutta type single step schemes.
The methods RKF43, RKF34, RKF54a, RKF54b, RKF45a, RKF45b, RKF87, RKF78, DOPRI54, DOPRI45, DOPRI65, DOPRI56, DOPRI87, DOPRI78, CK54, CK45 are embedded pairs of Runge-Kutta-Fehlberg, Dormand-Prince and Cash-Karp type, respectively. Estimates of the local discretization error are obtained in the usual way by comparing the results of the two submethods of the pair. The names indicate the orders of the subprocesses. For instance, RKF34 and RKF43 denote the same embedded Runge-Kutta-Fehlberg pair with orders 3 and 4. In RKF34 the result of the fourth order submethod is accepted, whereas RKF43 advances using the result of the third order submethod. In both cases the discretization error of the lower order subprocess is estimated and controlled.
For the single methods EULER1 (the first order Euler scheme), RK4 (the classical fourth order Runge-Kutta scheme) and BUTCHER6 (a Runge-Kutta scheme of order 6), the relative local error is controlled by comparing steps with different step sizes. The effective order of these methods is increased by one through local extrapolation.
Local extrapolation is also available for the submethods of the embedded pairs. For instance, the method xRKF78 uses extrapolation of the 8-th order subprocess of RKF78, yielding a method of effective order 9. The 7-th order subprocess is ignored. The cheap error estimate based on the embedded pair is not used implying some loss of efficiency when comparing RKF78 (order 8) and xRKF78 (effective order 9).
The call numeric::butcher(method) returns the Butcher data of the methods used in numeric::odesolve. Here method is one of EULER1, RKF43, RK4, RKF34, RKF54a, RKF54b, DOPRI54, CK54, RKF45a, RKF45b, DOPRI45, CK45, DOPRI65, DOPRI56, BUTCHER6, RKF87, DOPRI87, RKF78, DOPRI78.
Note: Only local errors are controlled by the adaptive mechanism. No control of the global error is provided! |
Note: The run time of the numerical integration with a method of order p grows like , when DIGITS is increased. Computations with high precision goals are very expensive! High order methods such as the default method DOPRI78 should be used. |
Presently, only single step methods of Runge-Kutta type are implemented. Stiff problems cannot be handled efficiently with explicit methods such as the default method DOPRI78. For stiff problems, one may use one of the implicit A-stable methods GAUSS(s). See Example 7.
For problems of the special type with a matrix valued function f(t, Y), there is a specialized ("geometric") integration routine numeric::odesolveGeometric. Generally, numeric::odesolve is faster than the geometric integrator. However, odesolveGeometric preserves certain invariants of the system more faithfully.
J.C. Butcher: "The Numerical Analysis of Ordinary Differential Equations", Wiley, Chichester (1987).
E. Hairer, S.P. Norsett and G. Wanner: "Solving Ordinary Differential Equations I", Springer, Berlin (1993).