Getting Consistent Initial Conditions for a DAE

Math à la Carte Report 2008-2

February, 2008 Fred T. Krogh and Richard J. Hanson

Abstract We propose here a method for getting consistent initial conditions for an index 1 Differential Algebraic Equation. This is to replace the current approach used in the code DDASLX which is a modification of Petzold’s DASSL.

Contents

1 Introduction 1

2 When Second Partials are Small 2 2.1 The Dense Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 The Sparse Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

3 When Second Partials Can Not be Ignored 3

1 Introduction

Given a differential algebraic equation of the form

f(t, y, y′) = 0  (1)

known values for some of the initial conditions may not be known. We assume that y is known and that some of the y′ are unknown. In order to get off to a good start these should be selected in a way that makes things consistent before starting the integration.

In thinking about this problem we have used the following example from [1, p. 138] in which the derivatives are unknown when starting.

\begin{aligned}y_{1} + y′_{1} + y′_{2} − g_{1}(t) & = 0 \\ y_{2} − g_{2}(t) & = 0 \end{aligned}  (2)

In this example it is clear that only the sum of the derivatives is determined by the first equation, and none of the derivatives appears in the second equation. Thus the usual Newton method for solving a system of equations will not work. It is also clear from the second equation that consistency would require that y′_{2} = g′_{2}(t).

The problem is to have an approach that works in the general case. In the index one case we require f_{t} which is not usually required by a DAE solver. This must be provided for in a code even though in many cases it will be 0.

2 When Second Partials are Small

We follow the approach given in [1, p. 138-139], but provide more detail on how to proceed. We want to satisfy the following two equations.

f(t, y, y′) = 0
\frac{d}{dt}f = f_{t} + f_{y}y′ + f_{y′}y′′ = 0

(3)

Using our example problem for guidance, we would expect this system to be rank deficient and under determined.

We assume we can ignore second partials of f with respect to y′ and thus solve for correction to y′ using the following system.

[\begin{matrix}f_{y′} & 0 \\ f_{y} & f_{y′} \end{matrix}] [\begin{matrix}δy′ \\ y′′ \end{matrix}] = [\begin{matrix}−f(t, y, y′) \\ −f_{t} − f_{y}y′ \end{matrix}]  (4)

2.1 The Dense Case

To solve this system we can use Gaussian elimination with full pivoting on f_{y′} until we no longer have a sufficiently large pivot. If f_{y′} is fully factored we can simply solve for y′ and f_{y} need not even be computed. Otherwise apply the same transformations in the bottom block as were applied in the top block. If in a given row we were not able to find a pivot in the top block, then the corresponding row in the bottom block we may assume to be 0 in the right bottom block. Thus if there is a solution there must be a pivot available in the left bottom block in this row, and this we could use. Otherwise the system is not of index 1. Note that this solution does not even require that we solve for y′′.

As a gross check on what all this means let’s apply this to our example problem. For this case our equations look like

[\begin{matrix}1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 \\ 0 & 1 & 0 & 0 \end{matrix}] [\begin{matrix}δy′_{1} \\ δy′_{2} \\ y′′_{1} \\ y′′_{2} \end{matrix}] = [\begin{matrix}g_{1} − y_{1} − y′_{1} − y′_{2} \\ g_{2} − y_{2} \\ g′_{1} − y′_{1} \\ g′_{2} − y′_{2} \end{matrix}]  (5)

Thus we would use the top row, and the bottom row of this matrix in the process of getting a solution for y′.

If f is linear in y′ no iteration is required. We know of no examples where this is not the case. If the nonlinear dependence on y′ is weak, then an iterations is likely to converge. Such an iteration would only need to do the factorization on the first iteration. Finally if such dependence is strong we offer Section 3 as a possible approach.

At the expense of getting a lot messier, the approach here could be extended to the case when some mixture of at most n of the y_{i} and y_{i}′ are unknown where n is the size of the system.

2.2 The Sparse Case

A black box sparse solver may have difficulty solving Eq. (4) as it stands. But since presumably the actual solution for y′′ is not important, we can add a condition that tends to minimize its norm. Imagine that we add a block of rows to this problem of the form εf_{y′}δy + εIy′′ ≅ 0. If we apply Given rotations to 0 out the left block from this newly added row using the top block of rows, the f_{y′} is multiplied by (1 + ε^{2})^{1/2}, and the 0 is replaced by ε^{2}/(1 + ε^{2})^{1/2}I. We then discard the bottom rows just added. and divide the top block of rows by (1 + ε^{2})^{1/2} the resulting system is the same as Eq. (4), except for the multiplication replacing the 0 by ε^{2}/(1 + ε^{2})I, and the original right hand side divided by (1 + ε^{2}). Of course if ε is slightly less than the square root of the machine epsilon, the only change is in the 0. This change should be enough to allow a direct sparse solver to get a solution.

3 When Second Partials Can Not be Ignored

If the above iteration does not converge it is probably because the second partial derivatives of f with respect to y′ are too large to ignore. One could compute numerical approximations to f_{yy′}y′, f_{y′y′}y′′, and f_{ty′}. But this requires more work and a lot more space that we would like to require. We suggest something else here which has some chance of working. Given that we know of no problems that require this, we have not tried it.

Our plan is to pick a direction and get the extra derivative information only in that direction. Let ỹ = Qy for some orthogonal matrix Q. Now we would like to include information that comes from the second partial with respect to ỹ_{1}. We start by getting the full partial with respect to ỹ.

\frac{∂}{∂y′} [f_{y}y′ + f_{y′}y′′ + f_{t}] = \frac{∂}{∂ỹ′} [f_{y}y′ + f_{y′}y′′ + f_{t}] \frac{∂ỹ′}{∂y′}

= \frac{∂}{∂ỹ′} [f_{y}Q^{T} ỹ′ + f_{y′}Q^{T} ỹ′′ + f_{t}] Q

= [\frac{∂}{∂ỹ′} f_{y}] Q^{T} ỹ′Q + f_{y} + [\frac{∂}{∂ỹ′} f_{y′}] Q^{T} ỹ′′Q + [\frac{∂}{∂ỹ′} f_{t}] Q

= [\frac{∂}{∂ỹ′} f_{y}] y′Q + f_{y} + [\frac{∂}{∂ỹ′} f_{y′}] y′′Q + [\frac{∂}{∂ỹ′} f_{t}] Q

(6)

The partial with respect to y′′ of course has not changed, but we have three new terms for the partial with respect to y′. With any luck we may discover that (∂/∂ỹ′)f_{t} is 0 and can be ignored. Otherwise this one will require forming differences of differences.

At this point we wish to consider only the contribution to these additional terms with respect to ỹ_{1}. We can compute approximations to ∂f_{y}/∂ỹ′_{1}, ∂f_{y′}/∂ỹ′_{1}, and (∂/∂ỹ′_{1})f_{t} by perturbing y′ a small amount in the direction of ỹ_{1}, which in the original coordinate system will be some direction d which we select. Denote the first two difference quotients i.e. the difference the partials divided by the magnitude of the difference in d by d_{y} and d_{y′}, which are full matrices. The products d_{y}y′ and d_{y′}y′′ are best thought of as full matrices which are 0 in every column but the first. Thus these matrices can be multiplied on the right by Q, giving us a better approximation to the full derivative with respect to y′. The contribution from the f_{t} term should also be thought of as a matrix with all columns 0 except for the first and thus this matrix can also multiply Q.

For the direction d we propose to use the change in y′ from the previous iteration, and select Q so that Qd is 0 in every component but the first. Another possible choice would be to use y′′ as this defines the direction in which we are planning to change y′ when taking a step.

References

[1] Campbell S. L. Brenan K. E. and Petzold L. R., Numerical Solution of Initial Value Problems in Differential Algebraic Equations, SIAM, Philadelphia (1996) 256 pages.