numeric::odesolve -- numerical
solution of an ordinary differential equation
Introductionnumeric::odesolve(t0..t, f, Y0, ..)
returns a numerical approximation of the solution Y(t) of
the first order differential equation (dynamical system)
dY/dt = f(t,Y), Y(t0) = Y0, t0 and t real, Y and Y0 complex vectors.
Call(s)numeric::odesolve(t0..t, f, Y0 <, method> <, RelativeError = tol> <, Stepsize = h> <, Alldata =
n> <, Symbolic> )
Parameterst0 |
- | numerical real value for the initial time |
t |
- | numerical real value (the ``time'') |
f |
- | procedure representing the vector field |
Y0 |
- | list or 1-dimensional array of numerical values representing the initial condition Y0 |
Optionsmethod |
- | name of the numerical scheme |
RelativeError =
tol |
- | forces internal numerical Runge-Kutta steps to use
stepsizes with local discretization errors below tol. This
tolerance must be a numerical real value >=
10^(-DIGITS). |
Stepsize = h |
- | switches off the internal error control and uses a
Runge-Kutta iteration with constant stepsize h.
h must be a numerical real value. |
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. |
Symbolic |
- | makes numeric::odesolve return a vector
of symbolic expressions representing a single symbolic step of the
Runge-Kutta iteration. |
ReturnsThe 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 Y0.
With the option Alldata a list of mesh data is
returned.
Side
EffectsThe function is sensitive to the environment variable DIGITS, which determines the
numerical working precision.
Related
Functionsnumeric::butcher,
numeric::odesolve2,
plot::ode
Detailsnumeric::odesolve is a general purpose solver able to
deal with initial value problems of various kinds of (non-stiff)
ordinary differential equations. Equations d^p/dt^p
y(t)=f(t,y,y',..) of order p can be solved by
numeric::odesolve after reformulation to dynamical system
form. This can alway be achieved by writing the equation as a first
order system
d/dt [Y[1],..,Y[p-1],Y[p]] = [Y[2],..,Y[p],f(t,Y[1],..,Y[p])]
for the vector [Y[1],Y[2],..]=[y,y',..]. Cf. example 4.t0, t and Y0
must not contain symbolic objects which cannot be converted to floating
point values via float.
Numerical expressions such as exp(PI),
sqrt(2) etc. are accepted.f defining the dynamical system
Y'=f(t,Y) 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
PI, exp(2) etc. However, all values must be
convertible to real or complex floating point numbers by float.t and Y.Y has to be represented by a
list or an array with one element. For instance, the input data for the
scalar initial value problem diff(y(t),t)=t*sin(y(t)),
y(0)=1 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 */ |
DIGITS: an adaptive
control of the stepsize keeps local relative discretization errors
below tol=10^(-DIGITS), unless a different
tolerance is specified via the option RelativeError=tol. The error control may be
switched off by specifying a fixed Stepsize=
h.Only local errors are controlled by the adaptive mechanism. No control of the global error is provided!
Y := t -> numeric::odesolve(t0..t, f, Y0
<,options>) the numerical solution can be repesented by a
MuPAD function: the call Y(t) will start the
numerical integration from t0 to t. A more
sophisticated form of this function may be generated via
Y := numeric::odesolve2(f, t0,
Y0 <, options>)
This equips Y with a remember mechanism that uses
previously computed values to speed up the computation. Cf.
example 2.
numeric::odesolveGeometric which preserves
geometric features of the system more faithfully than
numeric::odesolve.
Option:
method| 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), | BUTCHER6 (order 6), | RKF87 (order 7), |
| xRKF87 (order 7), | RKF78 (order 8), | xRKF78 (order 8). |
The Background below provides some information on these methods.
Option: RelativeError = toltol.tol = 10^(-DIGITS)
ensures that the local discretization errors are of the same order of
magnitude as numerical roundoff.
Usually there is no need to use this option to change this setting.
However, occasionally the numerical evaluation of the Runge-Kutta steps
may be ill-conditioned or stepsizes are so small that the time
parameter cannot be incremented by the stepsize within working
precision. In such a case this option may be used to bound the local
discretization error by tol and use a higher working
precision given by DIGITS.
tol >=
10^(-DIGITS) are accepted.The global error of the result returned by
numeric::odesolve is usually larger than the local errors
bounded by tol. Although the result is displayed with
DIGITS decimal places
one should not expect that all of them are correct. The relative
precision of the final result is tol at best!
Option: Stepsize = hnumeric::odesolve uses an adaptive
stepsize control mechanism in the numerical iteration. The option Stepsize=h switches off this adaptive
mechanism and uses the Runge-Kutta method specified by
method (or the default method RKF78)
with constant stepsize h.t of the integration interval t0..t, if
(t-t0)/h is not integer.With this option there is no automatic error control! Depending on the problem and on the order of the method the result may be a poor numerical approximation of the exact solution.
Option: Alldata = nnumeric::odesolve returns a list of
numerical mesh points
[[t0, Y0], [t1, Y1], .., [t, Y(t)]]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 [t0,Y0] and the final point
[t,Y(t)]. For n<=0 only the data
[[t0,Y0],[t,Y(t)]] are returned.
Option: Symbolicnumeric::odesolve(t0..t, f, Y0,
<method>, Symbolic) returns a vector
(list or array) of expressions representing a single step of the
numerical scheme with stepsize t-t0. In this mode symbolic
values for t0, t and the components of
Y0 are accepted.
Example
1We compute the numerical solution y(10) of the initial value problem y'=t*sin(y), y(0)=2:
>> f := proc(t, Y) begin [t*sin(Y[1])] end_proc:
>> numeric::odesolve(0..10, f, [2])
[3.141592654]
>> delete f:
Example
2We consider y'=y, y(0)=1. The numerical solution may be represented by the function
>> Y := t -> numeric::odesolve(0..t, (t,Y) -> Y, [1]):
Calling Y(t) starts the numerical
integration:
>> Y(-5), Y(0), Y(1), Y(PI)
[0.006737946999], [1.0], [2.718281828], [23.14069263]
>> delete Y:
Example
3We compute the numerical solution Y(PI)=[x(PI),y(PI)] of the system
x'=x+y, y'=x-y, x(0)=1, y(0)=I.
>> f := (t, Y) -> [Y[1] + Y[2], Y[1] - Y[2]]: Y0 := [1, I]:
>> numeric::odesolve(0..PI, f, Y0)
[72.57057162 + 30.05484302 I, 30.05484302 + 12.46088558 I]
The solution of a linear dynamical system Y'=A*Y with a constant matrix A is Y(t)=exp(t*A)*Y(0). 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)
[72.57057163 + 30.05484303 I, 30.05484303 + 12.46088558 I]
>> delete f, Y0, t, tA:
Example
4We compute the numerical solution y(1) of the differential equation y''=y^2 with initial conditions y(0)=0, y'(0)=1 The second order equation is converted to a first order system for Y=[y,y']=[y,z]:
y'=z, z'=y^2,y(0)=0, z(0)=1.
>> f := proc(t, Y) begin [Y[2], Y[1]^2] end_proc: Y0 := [0, 1]:
>> numeric::odesolve(0..1, f, Y0)
[1.087473533, 1.362851121]
>> delete f, Y0:
Example
5We demonstrate how numerical data can be obtained on a
user defined time mesh t[i]. The initial value problem is
y'=sin(t)-y, 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(t[i-1]..t[i], f, Y[i-1])
end_for:
The following mesh data are produced:
>> [t[i], Y[i]] $ i = 0..100
[[0, [1]], [1/10, [0.9097]], [1/5, [0.8374]], [3/10, [0.7813]],
[2/5, [0.7397]], ... , [99/10, [0.2159]], [10, [0.1476]]]
These data can be displayed by a list plot:
>> plotpoints := [point(t[i], op(Y[i])) $ i = 0..100]:
>> plot2d([Mode = List, plotpoints])
The same plot can be obtained directly via plot::ode:
>> plot(plot::ode(
[t[i] $ i = 0..100], f, Y[0],
[(t, Y) -> [t, Y[1]], Style = Points]))
>> delete f, t, DIGITS, Y, plotpoints:
Example
6We compute the numerical solution y(1) of y'=y, y(0)=1 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(0..1, f, [1], RK4)
[2.718281828459]
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(0..1, f, [1], xRKF78)
[2.718281828459]
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:
Example
7We consider the initial value problem y'=-10^(10)*y*(1-cos(y)), 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^10*Y[1]*(1 - cos(Y[1]))]: Y0 := [1]:
We first attempt to compute y(1) with a
working precision of 4 digits using the default setting
RelativeError =10^(-DIGITS)=10^(-4):
>> DIGITS := 4: numeric::odesolve(0..1, f, Y0)
[2.931e-5]
Due to numerical roundoff 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(0..1, f, Y0, RelativeError = 10^(-4))
[0.0000099999997577602193132]
Since relative local discretization errors are of the magnitude 10^(-4), 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(0..1, f, Y0, RelativeError = 10^(-10))
[0.000010000000000493745906]
>> delete f, Y0, DIGITS:
Example
8We compute the numerical solution y(1) of y'=y, y(0)=1 with various methods and various constant stepsizes. We compare the result with the exact solution y(1)=exp(1)=2.718281828....
>> f := (t, Y) -> Y: Y0 := [1]:
We first use the Euler method of order 1 with two different stepsizes:
>> Y := numeric::odesolve(0..1, f, Y0, EULER1, Stepsize = 0.1):
>> Y, globalerror = float(exp(1)) - Y[1]
[2.59374246], globalerror = 0.1245393684
Decreasing the stepsize by a factor of 10 should reduce the global error by a factor of 10. Indeed:
>> Y := numeric::odesolve(0..1, f, Y0, EULER1, Stepsize = 0.01):
>> Y, globalerror = float(exp(1)) - Y[1]
[2.70481383], globalerror = 0.01346799904
Next we use the classical Runge-Kutta method of order 4 with two different stepsizes:
>> Y := numeric::odesolve(0..1, f, Y0, RK4, Stepsize = 0.1):
>> Y, globalerror = float(exp(1)) - Y[1]
[2.718279744], globalerror = 0.000002084323879
Decreasing the stepsize 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(0..1, f, Y0, RK4, Stepsize = 0.01):
>> Y, globalerror = float(exp(1)) - Y[1]
[2.718281828], globalerror = 0.0000000002246438649
>> delete f, Y0, Y:
Example
9We integrate y'=y^2, y(0)=1 over the interval t=0..0.99 with a working precision of 4 digits. The exact solution is y(t)=1/(1-t). 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(0..0.99, f, Y0, Alldata)
[[0.0, [1.0]], [0.5668, [2.308]], [0.7784, [4.513]],
[0.8842, [8.636]], [0.9371, [15.9]], [0.9636, [27.43]],
[0.9768, [43.05]], [0.99, [99.99]]]
With Alldata=2, only
each second point is returned:
>> numeric::odesolve(0..0.99, f, Y0, Alldata = 2)
[[0.0, [1.0]], [0.7784, [4.513]], [0.9371, [15.9]],
[0.9768, [43.05]], [0.99, [99.99]]]
One can control the time mesh using the option Stepsize=h:
>> numeric::odesolve(0..0.99, f, Y0, Stepsize=0.1, Alldata = 1)
[[0.0, [1.0]], [0.1, [1.111]], [0.2, [1.25]], [0.3, [1.429]],
[0.4, [1.667]], [0.5, [2.0]], [0.6, [2.5]], [0.7, [3.333]],
[0.8, [5.0]], [0.9, [10.0]], [0.99, [131.2]]]
However, with the option Stepsize=h, no automatic error control is
provided by numeric::odesolve. Note the poor approximation
y(t)=131.1 for t=0.99 (the exact value is
y(0.99)=100.0). The next computation with smaller stepsize
yields better results:
>> numeric::odesolve(0..0.99, f, Y0, Stepsize = 0.01, Alldata = 10)
[[0.0, [1.0]], [0.1, [1.111]], [0.2, [1.25]], [0.3, [1.429]],
[0.4, [1.667]], [0.5, [2.0]], [0.6, [2.5]], [0.7, [3.333]],
[0.8, [5.0]], [0.9, [10.0]], [0.99, [100.0]]]
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:
Example
10The second order equation y''+sin(y)=0 is written as the dynamical system y'=z, z'=-sin(y) for the vector Y=[y,z]. A single symbolic step
[y(t0),z(t0)] -> [y(t0+h), z(t0+h)]
of the Euler method is computed:
>> f := proc(t, Y) begin [Y[2], -sin(Y[1])] end_proc:
>> numeric::odesolve(t0..t0+h, f, [y0, z0], EULER1, Symbolic)
[y0 + h z0, z0 - h sin(y0)]
>> delete f:
BackgroundJ.C. Butcher: ``The Numerical Analysis of Ordinary Differential Equations'', Wiley, Chichester (1987).
E. Hairer, S.P. Noersett and G. Wanner: ``Solving Ordinary Differential Equations I'', Springer, Berlin (1993).
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, BUTCHER6, RKF87, RKF78.Only local errors are controlled by the adaptive mechanism. No control of the global error is provided!
The run time of the numerical integration with a
method of order p grows like
O(10^(DIGITS/(p+1))), when DIGITS is increased. Computations
with high precision goals are very expensive! High order methods such
as the default method RKF78 should be used.
numeric::odesolveGeometric. Generally,
numeric::odesolve is faster than the geometric integrator.
However, odesolveGeometric preserves certain invariants of
the sytem more faithfully.tol
allows to set a precision goal independent of the working
precision.