|
Solution of the initial value problem of ode's
Adams-Bashforth-Moulton in PECE mode |
|
|
|
|
| |
- 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 grid 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 basic idea of this method is simple:
Using the new unknown value f(ti,yi) and a number of back
values f(ti-j,yi-j), j=1,..,k one constructs an interpolation polynomial
p(t) of order at most k which serves as an approximation for y'(t).
This interpolation polynomial then replaces y' in the formula
yi = yi-1 + ∫ti-1 to ti y'(t) d(t)
This results in an Adams-Moulton formula. Now in this formula the unknown value
yi appears as an argument of f, hence this is an implicit scheme.
The ''C'' in PECE stands for this (''corrector'').
Now one repeats the interpolation process, omitting the unknown value yi but
using the back values f(ti-j,yi-j), j=1,..,k again
and gets a second polynomial q of order k-1 (at most).
This again is inserted in the integral representation
of yi above, replacing y'. This defines a second approximation
for yi. The result is an explicit formula, Adams-Bashforth, and ''P'' (predictor) is the indicator for this.
This predictor value replaces the unknown yi inside f in the Moulton formula.
Finally f is reevaluated at the new point.
Hence this needs in all 2 function evaluations per integration step, for which the two ''E'' in PECE stand for.
The net result of this is a formula of the type
yi = yi-1 + ∑j=1 to k αi,j fi-j +
+αi,0f(ti,yi-1 +∑l=1 to k
βi,lfi-l)
fi = f(ti,yi)
with the αi,j the coefficients of the Moulton formula (i.e. the integrals of
the Lagrange polynomials over [ti-1,ti]) and βi,l
those of the Bashforth formula. These coefficients have two indices since the grid is completely
arbitrary. Rewriting the formulae using backward differences makes the recomputation of
the coefficients much easier.
- 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 theshhold is defined here as
reltol × ||y||(tcurrent)+abstol
with reltol, abstol userdefined.
One of the highlights of this approach is the fact that the difference of the predicted
and the corrected value for yi can serve as a estimator for this error
and that also the form of the local errors of both the predictor and the corrector is explicitly
known. The higher derivatives of y appearing in this formulae can be estimated safely
using divided differences of the f values.
- 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. The method sketched above has order k+1. Hence in principle you might
rise the order to any desired value using two function evaluations per step only.
- 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.
The method here is asymptotically stable for any choice of the grid and of the interpolation degree.
- 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 method here has a rather small region of absolute stability in the left half plane of varying shape,
containing a circle with center on the negative real axis and 0 as a boundary point with radius in the order
of 1, hence there must hold
hj | λ | <= C(k) with C(k) in the order of one.
However, C(k) shrinks as k increases, hence it makes no sense to use very large k
and this is limited to 12 here. This method is not adequate for stiff equations.
It is well suited for smooth solutions where high precision is required.
-
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 1) and takes the order which maximizes
the gridsize. 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.
- Methods of higher order are generally more efficient than those of low order.
- The computational code behind this is ''ODE'' of Shampine and Gordon, available
from NETLIB/ODE.
|
| |
|
|
| |
|
Input: Predefined cases |
| |
- There is a choice between 6 different problems.
- 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 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 local error is required to be
less or equal than (||y||+smally)*tol.
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 | | hmax | | tol | | smally | | num |
| 1 | | 0.5 | | 1.e-6 | | 0.1 | | 1 |
| 2 | | 0.2 | | 1.e-6 | | 1.0 | | 1 |
| 3 | | 0.1 | | 1.e-4 | | 1.e-6 | | 1 |
| 4 | | 0.25 | | 1.e-4 | | 0.1 | | 1 |
| 5 | | 0.1 | | 1.e-2 | | 1.e-8 | | 1 |
| 6 | | 0.3 | | 1.e-3 | | 1.0 | | 1 |
| 7 | | - | | - | | - | | - |
- For all cases h_min=8.8e-16 × tcurrent is the smallest time step allowed. If the stepsize control goes below this,
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 smally the stepsize selection and the precision?
Try testcase 3!
|
| |
|
|
| |
|
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 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 resulting in
reltol=tol and abstol=smally*tol .
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.
|
| |
|
|
| |