| |
- 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.
- All the methods we are using here restrict themselves 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 basic idea of these methods can be twofold: either one discretizes
y' at a time instance ti using neighboring already known
values and inserts this value in the differential equation, using the value
of the right hand side at ti too or one uses integration of the ode
over some time interval [ti-j,ti] and replaces
the right hand side integral by a quadrature rule using back (and possibly the unknown new)
values as values of f on the nodes.
Hence these methods are called ''integrators'' sometimes.
- ti denotes the new time point here, and there holds j=1,
but the number of back values used for interpolation varies. Since the methods used here
all result in an equation for the new grid value yih
appearing also inside f they are called ''implicit''. For nonlinear
f these equations are nonlinear and must be solved by iteration,
which of course is terminated after a finite number of steps. Hence strictly
speaking you can write this as an explicit formula for the new value, finally.
- 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 treshhold is defined here as
reltoli × |yi(tcurrent)|+abstoli, i=1,...,n
with reltol, abstol 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 methods which we present here are:
- The Adams-Moulton method of variable order. This is obtained in the following steps:
- step 1: rewrite the differential equation as
y(ti) = y(ti-1)+ ∫ti-1ti f(t,y(t)) dt
- step 2: replace in this integral f(t,y(t)) by a polynomial in t which interpolates back values including the
new one:
pk(tj) = f(tj,yhj) def= fj, j=i,i-1,..,i-k
Here yh denotes the computed solution. k varies from 1 to 11, yielding a consistency from
2 to 12.
- step 3: Integrate this polynomial exactly. Since the grid is not equidistant, this implies a recomputation
of the basis polynomials of the interpolation scheme. Therefore here the representation by backward differences is used:
pk(t) = ∑ j= 0 to k (
∏l=0 to j-1 (t-tl) ) ∇j fi
This results in a formula
yhi = yhi-1+ ∑j=0 to k αi,i-jfi-j.
- step 4: The equation for yhi is an implicit one, since this also accurs in fi.
This equation is now solved by direct iteration in the form
yhi,s = αi,i f(ti,yhi,s-1) + Gi ,
G denoting the term which is fixed in this iteration. The number m of iteration steps is limited to 4
in this code, otherwise a convergence failure is reported and the grid is refined.
The first guess is computed from an interpolation technique using only back values, the so called predictor (P).
Then any evaluation of f in the iteration counts as ''E'', and the computation of the next
yhi,s as a corrector (''C''). For the finally accepted value f
is evaluated again. In short this approach is named P(EC)mE, hence PECE for m=1.
Remark: for k=1 and exact solution of the implicit equation this would be the trapezoidal rule.
- This method is consistent for order k+1 for any k (we use no more than 12) and asymptotically stable, but has a rather small region of absolute
stability. On the real axis the interval is of the form [-C,0[ with C strongly decreasing with k increasing.
It is well suited for smooth solutions where high precision is required, but not useful for ''stiff'' equations
(for this see below).
- The BDF method, also with variable order and stepsize. It is obtained in the following manner:
- step 1: compute an interpolation polynomial of order <= k which satisfies
pk(ti-j) = yhi-j, j=0,..,k .
Observe that the value for j=0 is unknown.
- step 2: Differentiate this polynomial once and equate the value at ti with fi:
This results in an implicit equation
∑j=0 to k αi,i-j ∇j yhi = f(ti, yhi) .
- step 3: This equation is solved by the simplifed Newton's method (also known as ''chord iteration'') with
a Jacobian of f w.r.t. y, which is held piecewise constant (as long as convergence is
considered as sufficiently good).
- This method is consistent and asymptotically stable for k=1,...,5 with order k
and A-stable (the region of absolute stability covers the left complex halfplane) for k=1,2.
For k=3,4,5 it has a cone-like region of absolute stability in the left halfplane with the
opening angle decreasing with increasing k. It is therefore well suited for stiff equations.
An ode is said to be ''stiff'' if its solution manifold has components with extremely large growing
derivatives and simultaneously also very smooth components. Typically the smooth component is the
sought one, but the presence of the others hinders the use of adequate large stepsizes. Hence this is
is a term of a more qualitative than quantitative nature. A typical example:
y'(t) = -1000 (y(t) - exp(-t))
with the solution
y(t) = c exp(-1000 t)+(1000/999)exp(-t) .
The initial value fixes c but with t increasing any
solution tends to the smooth component (1000/999)exp(-t).
- In both cases the local error is analytically known with the form
constant × hp+1 × y(p+1)(ξ)
with ξ in the current integration/differentiation range and
p the order of consistency.
The higher derivative of the true solution
can be reliably estimated from divided differences of the back values of f.
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 in our output.
- 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.
- Discontinuities 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.
|