[Next Section] [Previous Section] [Table of Contents]

Full Equations (FEQ) Model for the Solution of the Full, Dynamic Equations of Motion for One-Dimensional Unsteady Flow in Open Channels and Through Control Structures

U.S. GEOLOGICAL SURVEY WATER-RESOURCES INVESTIGATIONS REPORT 96-4240

9.2 Solution of a Sparse, Banded Matrix


Repeated solution of the linear equation system by Newton's method results in a sequence of corrections that decrease to an acceptable value. Criteria are needed to determine when the corrections are small enough to end the iterations. A variety of problems may arise in the search for a solution and techniques for improving the convergence to a solution, so an evaluation of the solutions to these problems is needed.

The coefficient matrix for Newton's method for the system of nonlinear equations describing a stream network is banded, but the bandwidth is not constant. The equations for a branch have a constant bandwidth, but the algebraic relations for the special features add variability to the bandwidth at every junction (see section 10.3). The computation time for the solution of the linear equation system produced with Newton's method is reduced to a feasible value in FEQ by taking advantage of the band structure. A direct method based on the Crout variant of Gaussian elimination is used in FEQ simulation. The method as outlined here is derived from Zienkiewicz and Taylor (1989, p. 479-483).

The system of equations to be solved is

(133)

Equation ,

where

J is a Jacobian matrix of coefficients for the linear system,
u is the vector of unknowns, and
b is the right-hand-side vector of residuals.

To solve this system efficiently, let J = LU, where L is a lower triangular matrix and U is an upper triangular matrix. A triangular matrix is a square matrix in which all elements above or below the main diagonal are zero. Requiring the values along the main diagonal of L to be unity makes the two triangular matrices unique. To solve the system of linear equations, let Uu = w. The unknowns are then found by first computing w by use of the process of forward elimination. Once w is known, the next step, back substitution, solves for u as follows:

Equation

(134)
Equation

Equation

(135)
Equation .

The factoring of the matrix can be done in place; that is, the elements of the coefficient matrix are replaced by the factor of L or U as computed because the diagonal elements of L and the parts of L and U known to be zero are not stored. In-place factoring is assumed in the following equations. Thus, the same subscripts on any of J, L, or U refer to the same physical location in the stream system. The value that is needed will be at the location when reached. A row of L to the left of the diagonal (the diagonal element is 1.0 by definition) is calculated in the computations, followed by the column of U at and above the main diagonal. This is done for each row of L and column of U until all are computed; that is, for Equation ,

Equation

(136)

Equation .

The best way to visualize these operations is to partition the coefficient matrix into three parts: the reduced zone where L and U are already computed, the active zone that is currently computed, and the unreduced zone of the coefficient matrix. These zones are shown in figure 19.

These equations are written with the assumption that all the equations are stored with all the zero and nonzero coefficients. However, the vector inner product for the element of L in the j th row and the i th column involves only elements to the left of the i th column in the j th row and in the i th column at and above the main diagonal, as shown in figure 20. The pattern for the vector inner product for the element of U in the i th row of the j th column includes only elements to the right of the main diagonal in the i th row and above the i th row in the j th column, as shown in figure 21. The banded structure of the Jacobian matrix as developed previously is preserved with this pattern. This preservation of the banded structure results only if no row or column interchanges are made to prevent a zero divide in equation 136. The zeros within the band limits may become nonzero, but the zeros outside the band limits must remain as zeros. Therefore, both zero and nonzero coefficients, within the band as shown previously, are stored.

The band of the Jacobian matrix has constant width and pattern within branches. The junctions cause the band of the Jacobian matrix to be of variable width (section 10.3). Therefore, if enough nodes are in a branch, explicit computation of the vector inner products without iteration is worthwhile; that is, the pattern of the summation is constant for each element, so the summation is written explicitly. This technique is known as loop unrolling in algorithm design. Two block types are available in FEQ computations: a variable-band block and a fixed-band block for the branches. Use of loop unrolling greatly reduces solution time for factoring the matrix so that a fixed-band matrix can be factored with the variable-band algorithm at nearly the same speed as for a fixed-band algorithm.


[Next Section] [Previous Section] [Table of Contents]