|
Solution of the initial value problem of ode's
Gragg-Bulirsch-Stoer |
|
|
|
|
| |
- Any explicit ordinary differential equation, also such including higher order derivatives,
can be rewritten in the form
y'(t) = f(t,y(t))
Here y und f are elements in Rn, mostly n > 1.
y'(t) denotes the componentwise derivative of y(t) w.r.t
''time'' t. Hence we will use this ''standard'' form here.
One also may remove the special role of t by introducing a component
yn+1 by addition of the following equations to the
initial problem
y'n+1 = 1 ; yn+1(t0) = 0
substituting yn+1 for t, in which case the system is
called ''autonomous''.
- In order to get an uniquely solvable problem one needs n initial
values for y at some time point a=t0.
This means we fix some element from the solution manifold of this equation.
We may write this as
y(a) = ya.
In addition one needs beyond the joint continuity of f w.r.t. t,y some additional
properties, for example Lipschitz continuity w.r.t. y, at least locally
near (a,y(a)), in order to guarantee existence and uniqueness of the solution.
- The method we are using here restricts itself in finding approximate
values for y on a grid of the time axis. In the present case this is variable in
gridsize. The endpoint b=tend of this grid is fixed in advance, and of course one
assumes that the solution will exist over [a,b]. This is not trivial since
there are equations with singularities whose position depends also on the
initial value (for example the predefined testcase 2)
- The method presented here is based on the so called explicit midpoint rule and its
analytical properties for highly smooth functions f.
The explicit midpoint rule results form discretizing the equation
yi+1 = yi-1 +∫ti-1ti+1 f(t,y(t)) dt
using the simplest open Newton-Cotes formula for the integral
yhi+1 = yhi-1 + (ti+1-ti-1)*f(ti,yhi)
resulting in an explicit linear 2-step method. Here ti+1 denotes the next (new)
gridpoint and for initializing one needs not only the initial point but also some value
for gridpoint 1. usually this is obtained from an explicit Euler step.
We here use a modification of Kiehl and Zenger, which uses a (second order)
Taylor step for this, with the second derivative obtained by a finite difference formula
with a fine grid. The effect of this is discussed below. The upper index h discerns
the discretized from the true solution.
- Since the use of a predefined fixed grid size is often inefficient or even unsafe,
due to a possible blow up of the discretization error or an unnecessary large amount of function
evaluations, there is the possibility
to use a variably sized grid whose grid size is adapted locally, based on
additional computations. This is done by estimating the currently introduced
error (the so called ''local discretization error'') and one tries to keep
this one below an user defined threshhold by adapting the local gridsize.
This threshhold is defined here as a weighted Euclidean norm, which uses a componentwise relative error
replacing values below smally by smally, see below and with
smally userdefined.
- There are three different properties of such a method which are of importance:
consistency, asymptotic stability and absolute stability.
- One says the order of consistency of such a discretization is p if the
error of one step can be bounded from above by a constant times grid size p+1
locally.
- One says that a method is asymptotically stable if the computed solutions for two such
ode systems with neighboring functions f and g on a compact interval [a,a+δ]
say yf and yg, differ by no more than
exp(C*(δ)) × maxt in [a,a+δ] ||f(t,yf(t))-g(t,yg(t))||,
in words: a perturbation of the differential equation results in a perturbation of
the discrete approximate solution which is bounded by the size of this perturbation times a factor
growing exponential in time, for arbitrary small grid sizes.
Warning: there are many methods which are consistent to high order but not
asymptotically stable.
- One says that a method has a region G of absolute stability in the
complex plane, if the discretized solution of
y' = λ*y , y(0)=1 with Re(λ) < 0
using the fixed step size h decays in absolute value provided the
product h*λ is in G.
This property is important to allow rather large grid sizes for fast decaying stable systems if
one is not in the fast transient part of the solution.
- The explicit midpoint rule has order two, is asymptotically stable but has no interval
of absolute stability on the real axis, (but an interval of neutral stability from
-sqrt(-1) to sqrt(-1) on the imaginary axis). Hence it is as such almost
never used. It has however for locally fixed stepsize a remarkable property:
for sufficiently smooth f and in case of a constant stepsize h one can write
yhN = y(t) + ∑i=1 to m h2i *(
vi(t)+(-1)Nui(t) ) + O(h2m+2) (*)
where t-t0=Nh with smooth functions ui, vi
which however might grow with time t.
From this it becomes clear that the discretized solution might exhibit a growing oscillating
component. The use of the second order Taylor step as a first step eliminates u1.
(See Kiehl and Zenger in Computing 41, (1989), pages 131-136).
In addition one uses only even numbered grid points as approximate solution points in order
to avoid the nonsmooth oscillations. Then the form of the discretization error
suggests to get higher (even) order estimates of the true value y(t) by
the same technique as is used in Rombergs scheme for quadrature: at
t=ti we use a step
size sequence
h(j)=(ti+1-ti)/n(j) (here with n(j)=2j, j=1,2,3,..)
and compute the values
yh(j)i+n(j), j=1,...,m+1
all as approximations for y(ti+1).
Because of (*) we may consider these values as values of a (slightly perturbed)
polynomial in the variable h2 of degree m.
The value of this polynomial, evaluated at the argument h=0 produces its
absolut coefficient, hence y(ti+1), up to an error of size O(h(1)2m+2).
This evaluation is directly done using Neville's interpolation scheme (applied componentwise
for the vector y).
This gives a lower triangular array of approximations for y(ti+1) of order
2k in column k=1,2,...,m and in each row the difference of an element to its right neighbor
gives an (asymptotically correct) estimate for its (local) error: we get a cheap method of error
estimation as well. This technique is applied piecewise, on each subinterval of the grid.
- The local error bound with the estimation of the local error (which of course also depends on the current order used)
dictates the local gridsize. The code tries to vary the order (but by no more than 2) and takes the order which minimizes
an estimate of the total effort. A plot of the current order over time is displayed.
- Whereas the control of the local error works reasonably for the variable stepsize mode
there is no control for the global error which is determined by the properties of
the solution manifold of the ode.
- One observes strong variations of the stepsize if the derivatives of y (resp. f)
vary strongly.
- Discontinuites of f or of its higher partial derivatives
reduce the stepsize often dramatically. Hence it is advisable to identify their
locations and to restart the integration from there.
- ''Stiff'' equations (i.e. something corresponding to an huge negative
λ in the model equation above) let the explicit integrators fail.
Moreover, this method can produce instabilities if the stepsize is not
small enough, this being not necessarily detected by the integrator.
(Take testcase 4 with low and high precision requirement).
- This method restricts not the stepsize from below but limits
the number of successive stepsize reductions from the same time point
(there may be of course several) to 9. In case of failure you
get a corresponding message. This means that the equation is extremely
stiff. You may try to reduce the parameter hmax, but we limit
the number of steps to be taken to 100000 anyway.
- Methods of higher order are generally more efficient than those of low order,
but since this method tries to minimize the total effort you may observe that
in simple cases the order stays at its lower bound 2.
|
| |
|
|
| |
|
Input: Predefined cases |
| |
- There is a choice between 6 different problems which can be solved with the two
different integrators.
- These problems are
- the onedimensional model equation y'(t) = -2y(t)
with initial value y(0) = 1.
- the onedimensional test case y'(t) = (y(t))2
with initial value y(0) = 0.25.
- The restricted three body problem .
- computation of the current in a nonlinear RLC-circuit
with tunnel diode in parallel (simulating a Flip-Flop. A tunnel diode can act as a resistor or as a source,
depending on voltage).
- the reduced Reynolds equation.
- an ode with discontinuous right hand side
f(t,y(t)).
- You have the choice to limit the maximal attainable order between 2 and 16.
- You specify the maximal allowed stepsize to be taken (which of course must be considerably smaller than
tend-t0, t0=a above).
- In order to specify your bound for the local discretization error you specify
two parameters tol and smally.
The measurement of the local error is done in a weighted euclidean norm
||W*err||2 <= tol
where the components of the diagonal matrix W are defined as
1/max{|yhi|, |yhi-1|, smally}
and err is one of the difference vector in the Neville scheme.
Important: 0 < tol < 1 and smally > 0 !
- The parameter smally is a lower bound for the absolute value of the components of y
which play a role in the definition of the relative error.
- If you want the solution on a fixed predefined grid which might have no relation to the
internally used one you have the possibility to require output on an equidistant grid, defined
by an integer NUM (i.e. this grid size is (tend-t0)/NUM. It is not the integrators stepsize!)
- For the test cases 1,2 and 6 the exact solution is known, therefore we protocol the global error in these
cases.
- there are predefined values for the different cases as follows
| case | | hstart | | hmax | | tol | | smally | | num |
| 1 | | 0.1 | | 0.5 | | 1.e-6 | | 0.1 | | 1 |
| 2 | | 0.1 | | 0.2 | | 1.e-6 | | 1.0 | | 1 |
| 3 | | 0.001 | | 0.1 | | 1.e-4 | | 1.e-6 | | 1 |
| 4 | | 0.001 | | 0.25 | | 1.e-4 | | 0.1 | | 1 |
| 5 | | 0.000001 | | 0.1 | | 1.e-2 | | 1.e-8 | | 1 |
| 6 | | 0.03 | | 0.3 | | 1.e-3 | | 1.0 | | 1 |
| 7 | | 0.001 | | - | | - | | - | | - |
- For all cases the maximum number of integration steps is limited to 100000. If this is exceeded
then the integration is terminated, but you obtain the results computed so far.
|
| |
|
Output: predefined case |
| |
- The total computation time and an account of the effort used :
- the number of computed time steps.
- the number of rejected time steps (local error too large, step size reduced).
- the number of evaluations of f.
- In addition you get a plot of:
- step sizes used (in log form).
- local error estimation(in log form).
- the actual order used
- For the test cases 1,2 and 6 the exact solution is known,
therefore we protocol the global error in these cases.
(in log form of course)
- one or two components of the solution vector.
- !!!! It may occur that a log-scale displays no values
or a message like "all points undefined" appears.
this means that the errors are so small that by the default values of
the graphic program used these are all considered as zero.
This might occur with very strong precision requirements.
- Should the integrator issue warning messages you get these as a protocol.
|
| |
|
Questions: Predefined case?! |
| |
- Explain the behaviour of the stepsize for the different test cases.
- Behaves the local error as expected?
- In which cases is the global error indeed at most the sum the local errors?
- How influences the choice of tol resp. smally the stepsize selection and the precision?
Try testcases 3 and 4!
|
| |
|
|
| |
|
Input: user defined case |
| |
- Here you have the possibility to define an IVP yourself. This requires
- Specification of n. n is the number of components of
y resp. f(t,y).
- You write a piece of code for computing the components of f.
Details are given on the input form and in the help text.
- The initial values y(a).
- For graphical output you can specify two components for plotting.
- You have the choice to limit the maximal attainable order in the range 2 to 16 (even of course!)
- You specify the maximal allowed stepsize to be taken (which of course must be considerably smaller than
tend-t0, t0=a above.
- In order to specify your bound for the local discretization error you specify
two parameters tol and smally.
The measurement of the local error is done in a weighted euclidean norm
||W*err||2 <= tol
where the components of the diagonal matrix W are defined as
1/max{|yhi|, |yhi-1|, smally}
and err is one of the difference vector in the Neville scheme.
Important: 0 < tol < 1 and smally > 0 !
- The parameter smally is a lower bound for the absolute value of the components of y
which play a role in the definition of the relative error.
- If you want the solution on a fixed predefined grid which might have no relation to the
internally used one, you have the possibility to require output on an equidistant grid, defined by
by an integer NUM (i.e. this grid size is (tend-t0)/NUM. It is not the integrators stepsize!)
|
| |
|
Output: user defined case |
| |
- The total computation time and an account of the effort used :
- the number of computed time steps.
- the number of rejected time steps (local error too large, step size reduced).
- the number of evaluations of f.
- In addition you get a plot of:
- step sizes used (in log form).
- local error estimation(in log form).
- the order actually used
- one or two components of the solution vector.
- !!!! It may occur that a log-scale displays no values
or a message like "all points undefined" appears.
this means that the errors are so small that by the default values of
the graphic program used these are all considered as zero.
This occurs if the precision requirement is very strong, for example.
|
| |
|
Questions: user defined case?! |
| |
- This case is meant for the solution of an initial value problem of your
own interest.
- The questions to be asked here are of course the same as for the predefined cases.
|
| |
|
|
| |