numeric::odesolve
-- numerical
solution of an ordinary differential equationnumeric::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.
numeric::odesolve(t0..t, f, Y0 <, method> <, RelativeError = tol> <, Stepsize = h> <, Alldata =
n> <, Symbolic> )
t0 |
- | 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 |
method |
- | 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. |
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 Y0
.
With the option Alldata a list of mesh data is
returned.
The function is sensitive to the environment variable DIGITS
, which determines the
numerical working precision.
numeric::butcher
,
numeric::odesolve2
,
plot::ode
numeric::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
.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.
tol
.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!
numeric::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.
numeric::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.numeric::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.We 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:
We 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:
We 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:
We 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:
We 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:
We 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:
We 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:
We 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:
We 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:
The 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:
J.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.