|
Solution of the initial value problem for implicit ode's and dae's
Multistep method BDF (DASSL) |
|
|
|
|
| |
- In the following, y, f and G denote functions with
n components and y' denotes componentwise differentiation.
Typically where will be a ''free'' real parameter t with the
meaning ''time''. We search for a solution y as function of t,
given some relation between y, y' and t and some values of these
at some timeinstance t0. This is the field of
ordinary differential equations.
- In the applications an ordinary differential equation typically
does not appear in the explicit standard form y' = f(t,y) used by textbooks,
but in more general forms like
A(y)y' - f(t,y) = 0
or even more general
G(t,y,y') = 0 .
Here y' is defined only implicitly. Even more involved are cases
where not all components of the vector y' appear in these equations.
In this case the system is called a differential algebraic equation (DAE).
A famous example (and hard test case) is the constrained motion of a
pendulum:
x'' = λ x
y'' = λ y - g
0 = x2 + y2 - L
where λ (t) denotes the force acting along the pendulum and
x(t), y(t) denotes the position of its tip.
- In the following we assume that G is differentiable with
respect to all its variables at least twice.
-
There is an important notion in connection with DAE's, their index.
If the rank of the Jacobian of G with respect to y' is n
(the simple regular case) then the index is 0 and the system can be
solved pointwise, given t and y for y' using
e.g. Newtons method. In all other cases, especially if some components of
y' do not appear at all (as above for λ )
then these components are named algebraic components and the
index is 1 or greater. There are different notions of this ''index''.
One might try to transform the system by further differentiation finally
into a system of index 0. The minimum number of differentiation steps
required for this is the ''differential index'' used here.
The precise notion of this as follows: Consider the equation
G(t,y,y') = 0
and its total derivatives with respect to t:
(d/dt)j G(t,y,y') = 0 , j=1,...,k-1
and write this as a system of n × k equations
F(t,y,z) = 0 with z=(y',...,y(k)).
The smallest number k for which this equation determines
y' uniquely as a function of t and y is the
differential index.
Of course one gets higher order systems first in this way which however
finally might be transformed in a larger system of order one. There
remains finally the trouble to find initial values for the new solution
components so introduced, which is not an easy task. Dependent on the
rank of the joined Jacobian of G with respect to y and
y' it may be possible to solve the system
G(t0,y0,y'0) = 0
for the missing values. The code used here (DASKR from Petzold,
Brown, Hindmarsh and Ulrich) tries this, but for your tests you better did provide consistent
values yourself.
In practice only parts of the system need to be differentiated and by proper use
of the relations it might be possible to avoid additional dimensions.
- One must observe that for a proper DAE one doesn't have n
degrees of freedom in choosing the initial value. For example for the
pendulum example, there are only two degrees of freedom. If for example the system has
the semiexplicit form
F1(t,y,z,y') = 0 ; F2(t,y,z) = 0
with the Jacobian of F2 with respect to z of
full rank, then there remain clearly dim(F1) degrees of freedom only
in doing that.
- The basic idea of the method used here is the following:
The solution is constructed on a grid of the t-axis with
variable gridsize. Given the values of y at t0,...,
tm-1 one expresses y'(tm) using
the (unknown) value of y(tm) and back values
y(tm-1),..,y(tm-k) by a formula (''backward difference formula'')
∑i=0k αi,mym-i = hm-1 βk y'(tm)
inserts this into G(tm,ym,y'm)=0 and tries to solve now for ym.
The solver used for this is the simplified Newton's method (with only one Jacobian evaluation) and as
linear solver Gauss LU-decomposition with partial pivoting.
In the beginning one must use k=1 of course (this is known as ''Euler implicit'')
but from step to step k might be increased. This is the so called consistency
order of the formula, and for reasons of stability we use k <= 5 only.
(The case k=6 is theoretically possible, but has poor stability properties.)
Consistency order means that the error in y' is of the form
const × hk , h = max { hm,..,hm-k }
if back values would be exact.
- The local error of this method is controlled by computing the (analytically known)
error formula with the given data and stepsize and order are selected to get as large steps
as possible under the restriction to hold the local error below some user defined threshhold.
This theshhold 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 G and H on a compact interval [a,a+δ]
say yG and yH, differ by no more than
exp(C*(δ)) × maxt in [a,a+δ] ||G(t,yG(t),y'G)-H(t,yH(t),y'H(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.
- 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.
- 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.
- The convergence properties of these BDF methods with variable stepsize are only partially known,
for example one knows that for a linear system the order of convergence is not touched if
successive stepsize quotients are uniformly bounded.
- For DAE with index larger than one a new complication arises: here the algebraic components may
show up a convergence order less than the differential components, indeed the order may drop down to
k-index+2. This is known as ''order reduction''.
Therefore this application allows to exclude the algebraic components from the error control.
|
| |
|
|
| |
|
Input: Predefined cases |
| |
- There is a choice between 3 different problems.
- These problems are
- An externally controlled chemical reaction
- the van der Pol equation and its limit case
- the pendulum .
You have the possibility to use internally defined consistent initial values or to choose
free parameters (explained in the problem description) of the problem from which consistent
initial values are computed internally.
- You have the choice to limit the maximal attainable order, e.g. 2 in order to use A-stable
integrators only.
- 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
reltoli=tol and abstoli=smally*tol for all i
Important: 0 < tol < 1 and smally > 0 !
- The parameter smally resp. abstoli 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!)
- 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.0 | | 1 |
- For all cases h_min=1.0e-12 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 rejected steps because of slow convergence in the equation solving.
- the number of lu-decompositions.
- the number of evaluations of G.
- the number of evaluations of the joint Jacobian (dG(t,y,y')/d(y,y')).
- In addition you get a plot of:
- step sizes used (in log form).
- local error estimation(in log form).
- the actual order 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 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 omission of the algebraic components in the error test
effort and final precision?
|
| |
|
|
| |
|
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. G(t,y,y').
- You write a piece of code for computing the components of G.
Details are given on the input form and in the help text.
- The initial values y(t0), y'(t0).
- For graphical output you can specify two components for plotting.
- You have the choice to limit the maximal attainable order, e.g. 2 for BDF in order to use A-stable
integrators only.
- You specify the maximal allowed stepsize to be taken (which of course must be considerably smaller than
tend-t0, t0=a above.
- You also specify the initial used stepsize
- In order to specify your bound for the local discretization error you either specify
two parameters tol and smally resulting in
reltoli=tol and abstoli=smally*tol for all i
or you specify these parameters individually for each component of y.
Important: 0 < tol < 1 and smally > 0 !
- The parameter smally resp. abstoli 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 rejected steps because of slow convergence in the equation solving.
- the number of evaluations of G.
- the number of evaluations of the joint Jacobian (dG(t,y,y')/d(y,y')).
- the number of lu-decompositions.
- 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.
|
| |
|
|
| |