김경종
542 lines
28 KiB
Markdown
542 lines
28 KiB
Markdown
<!-- source-page: 941 -->
|
|
|
|
It should be noted that the Jacobi solution methods solve simultaneously for all eigenvalues and corresponding eigenvectors. However, in finite element analysis we require in most cases only some eigenpairs, and the use of a Jacobi solution procedure can be very inefficient, in particular when the order of K and M is large. In such cases much more effective solution methods are available, which solve only for the specific eigenvalues and eigenvectors that are actually required. However, the generalized Jacobi solution method presented in this section can be used very effectively as part of those solution strategies (see Section 11.6). When the order of the matrices K and M is relatively small, the solution of the eigenproblem is not very expensive, and the Jacobi iteration may also be attractive because of its simplicity and elegance of solution.
|
|
|
|
Subroutine JACOBI. Program JACOBI is used to calculate all eigenvalues and corresponding eigenvectors of the generalized eigenproblem $K\phi = \lambda M\phi$ . The argument variables and use of the subroutine are defined using comment lines in the program.
|
|
|
|
```csv
|
|
SUBROUTINE JACOBI (A,B,X,EIGV,D,N,RTOL,NSMAX,IFPR,IOUT) JAC00001
|
|
C ......
|
|
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . JAC00002
|
|
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . J AC00003
|
|
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . JAC 00004
|
|
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
A(N,N) = STIFFNESS MATRIX (ASSUMED POSITIVE DEFINITE) JAC00009
|
|
B(N,N) = MASS MATRIX (ASSUMED POSITIVE DEFINITE) JAC00010
|
|
X(N,N) = STORAGE FOR EIGENVECTORS JAC00011
|
|
EIGV(N) = STORAGE FOR EIGENVALUES JAC00012
|
|
D(N) = WORKING VECTOR JAC00013
|
|
N = ORDER OF MATRICES A AND B JAC00014
|
|
RTOL = CONVERGENCE TOLERANCE (USUALLY SET TO 10.**-12) JAC00015
|
|
NSMAX = MAXIMUM NUMBER OF SWEEPS ALLOWED JAC00016
|
|
(USUALLY SET TO 15) JAC00017
|
|
IFPR = FLAG FOR PRINTING DURING ITERATION JAC00018
|
|
EQ.0 NO PRINTING JAC00019
|
|
EQ.1 INTERMEDIATE RESULTS ARE PRINTED JAC00020
|
|
IOUT = UNIT NUMBER USED FOR OUTPUT JAC00021
|
|
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C .
|
|
A(N,N) = DIAGONALIZED STIFFNESS MATRIX JAC00024
|
|
B(N,N) = DIGONALIZED MASS MATRIX JAC00025
|
|
X(N,N) = EIGENVECTORS STORED COLUMNWISE JAC00026
|
|
EIGV(N) = EIGENVALUES JAC00027
|
|
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C .
|
|
IMPLICIT DOUBLE PRECISION (A-H,O-Z) JAC00030
|
|
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C . THIS PROGRAM IS USED IN SINGLE PRECISION ARITHMETIC ON CRAY JAC00032
|
|
EQUIPMENT AND DOUBLE PRECISION ARITHMETIC ON IBM MACHINES, JAC00033
|
|
ENGINEERING WORKSTATIONS AND PCS. DEACTIVATE ABOVE LINE FOR JAC00034
|
|
SINGLE PRECISION ARITHMETIC. JAC00035
|
|
DIMENSION A(N,N),B(N,N),X(N,N),EIGV(N),D(N) JAC00036
|
|
JAC00037
|
|
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 INITIALIZE EIGENVALUE AND EIGENVECTOR MATRICES JAC00038
|
|
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
```
|
|
|
|
<!-- source-page: 942 -->
|
|
|
|
```csv
|
|
20 X(I,J)=0. JAC00049
|
|
30 X(I,I)=1. JAC00050
|
|
IF (N.EQ.1) GO TO 900 JAC00051
|
|
C JAC00052
|
|
C INITIALIZE SWEEP COUNTER AND BEGIN ITERATION JAC00053
|
|
C JAC00054
|
|
NSWEEP=0 JAC00055
|
|
NR=N - 1 JAC00056
|
|
40 NSWEEP=NSWEEP + 1 JAC00057
|
|
IF (IFPR.EQ.1) WRITE (IOUT,2000) NSWEEP JAC00058
|
|
C JAC00059
|
|
C CHECK IF PRESENT OFF-DIAGONAL ELEMENT IS LARGE ENOUGH TO JAC00060
|
|
C REQUIRE ZEROING JAC00061
|
|
C JAC00062
|
|
EPS=(.01)**(NSWEEP*2) JAC00063
|
|
DO 210 J=1,NR JAC00064
|
|
JJ=J + 1 JAC00065
|
|
DO 210 K=JJ,N JAC00066
|
|
EPTOLA=(A(J,K)/A(J,J))*(A(J,K)/A(K,K))
|
|
EPTOLB=(B(J,K)/B(J,J))*(B(J,K)/B(K,K))
|
|
IF (EPTOLA.LT.EPS .AND. EPTOLB.LT.EPS) GO TO 210 JAC00069
|
|
C JAC00070
|
|
C IF ZEROING IS REQUIRED, CALCULATE THE ROTATION MATRIX JAC00071
|
|
C ELEMENTS CA AND CG JAC00072
|
|
C JAC00073
|
|
AKK=A(K,K)*B(J,K) - B(K,K)*A(J,K) JAC00074
|
|
AJJ=A(J,J)*B(J,K) - B(J,J)*A(J,K) JAC00075
|
|
AB=A(J,J)*B(K,K) - A(K,K)*B(J,J) JAC00076
|
|
SCALE=A(K,K)*B(K,K) JAC00077
|
|
ABCH=AB/SCALE JAC00078
|
|
AKKCH=AKK/SCALE JAC00079
|
|
AJJCH=AJJ/SCALE JAC00080
|
|
CHECK=(ABCH*ABCH + 4.*AKKCH*AJJCH)/4. JAC00081
|
|
IF (CHECK) 50,60,60 JAC00082
|
|
50 WRITE (IOUT,2020) JAC00083
|
|
GO TO 800 JAC00084
|
|
60 SQCH=SCALE*SQRT(CHECK) JAC00085
|
|
D1=AB/2. + SQCH JAC00086
|
|
D2=AB/2. - SQCH JAC00087
|
|
DEN=D1 JAC00088
|
|
IF (ABS(D2).GT.ABS(D1)) DEN=D2 JAC00089
|
|
IF (DEN) 80,70,80 JAC00090
|
|
70 CA=0. JAC00091
|
|
CG=-A(J,K)/A(K,K) JAC00092
|
|
GO TO 90 JAC00093
|
|
80 CA=AKK/DEN JAC00094
|
|
CG=-AJJ/DEN JAC00095
|
|
C PERFORM THE GENERALIZED ROTATION TO ZERO ELEMENTS JAC00096
|
|
C JAC00097
|
|
C JAC00098
|
|
90 IF (N-2) 100,190,100 JAC00099
|
|
100 JP1=J + 1 JAC00100
|
|
JMI=J - 1 JAC00101
|
|
KP1=K + 1 JAC00102
|
|
KMI=K - 1 JAC00103
|
|
IF (JM1-1) 130,110,110 JAC00104
|
|
110 DO 120 I=1,JMI JAC00105
|
|
AJ=A(I,J) JAC00106
|
|
BJ=B(I,J) JAC00107
|
|
AK=A(I,K) JAC00108
|
|
BK=B(I,K) JAC00109
|
|
A(I,J)=AJ + CG*AK JAC00110
|
|
B(I,J)=BJ + CG*BK JAC00111
|
|
A(I,K)=AK + CA*AJ JAC00112
|
|
120 B(I,K)=BK + CA*BJ JAC00113
|
|
130 IF (KP1-N) 140,140,160 JAC00114
|
|
140 DO 150 I=KP1,N JAC00115
|
|
AJ=A(J,I) JAC00116
|
|
BJ=B(J,I) JAC00117
|
|
AK=A(K,I) JAC00118
|
|
BK=B(K,I) JAC00119
|
|
```
|
|
|
|
<!-- source-page: 943 -->
|
|
|
|
```csv
|
|
A(J,I)=AJ + CG*AK
|
|
B(J,I)=BJ + CG*BK
|
|
A(K,I)=AK + CA*AJ
|
|
150 B(K,I)=BK + CA*BJ
|
|
160 IF (JP1-KM1) 170,170,190
|
|
170 DO 180 I=JP1,KM1
|
|
AJ=A(J,I)
|
|
BJ=B(J,I)
|
|
AK=A(I,K)
|
|
BK=B(I,K)
|
|
A(J,I)=AJ + CG*AK
|
|
B(J,I)=BJ + CG*BK
|
|
A(I,K)=AK + CA*AJ
|
|
180 B(I,K)=BK + CA*BJ
|
|
190 AK=A(K,K)
|
|
BK=B(K,K)
|
|
A(K,K)=AK + 2.*CA*A(J,K) + CA*CA*A(J,J)
|
|
B(K,K)=BK + 2.*CA*B(J,K) + CA*CA*B(J,J)
|
|
A(J,J)=A(J,J) + 2.*CG*A(J,K) + CG*CG*AK
|
|
B(J,J)=B(J,J) + 2.*CG*B(J,K) + CG*CG*BK
|
|
A(J,K)=0.
|
|
B(J,K)=0.
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
C
|
|
```
|
|
|
|
<!-- source-page: 944 -->
|
|
|
|
```csv
|
|
BB=SQRT(B(J,J))
|
|
DO 270 K=1,N
|
|
270 X(K,J)=X(K,J)/BB
|
|
GO TO 900
|
|
C
|
|
C UPDATE D MATRIX AND START NEW SWEEP, IF ALLOWED
|
|
C
|
|
280 DO 290 I=1,N
|
|
290 D(I)=EIGV(I)
|
|
IF (NSWEEP.LT.NSMAX) GO TO 40
|
|
GO TO 255
|
|
C
|
|
800 STOP
|
|
900 RETURN
|
|
C
|
|
2000 FORMAT (//,' SWEEP NUMBER IN *JACOBI* = ',I8)
|
|
2010 FORMAT ( ' ',6E20.12)
|
|
2020 FORMAT (//,' *** ERROR *** SOLUTION STOP',/,
|
|
1 ' MATRICES NOT POSITIVE DEFINITE')
|
|
2030 FORMAT (//,' CURRENT EIGENVALUES IN *JACOBI* ARE',/)
|
|
END
|
|
JAC00190
|
|
JAC00191
|
|
JAC00192
|
|
JAC00193
|
|
JAC00194
|
|
JAC00195
|
|
JAC00196
|
|
JAC00197
|
|
JAC00198
|
|
JAC00199
|
|
JAC00200
|
|
JAC00201
|
|
JAC00202
|
|
JAC00203
|
|
JAC00204
|
|
JAC00205
|
|
JAC00206
|
|
JAC00207
|
|
JAC00208
|
|
JAC00209
|
|
JAC00210
|
|
```
|
|
|
|
# 11.3.3 The Householder-QR-Inverse Iteration Solution
|
|
|
|
Another most important transformation solution technique is the Householder-QR-inverse iteration (HQRI) method, although this procedure is restricted to the solution of the standard eigenproblem (see J. G. F. Francis [A], J. H. Wilkinson [B], B. N. Parlett [A, B], and R. S. Martin, C. Reinsch, and J. H. Wilkinson [A]). Therefore, if the generalized eigenproblem $K\phi = \lambda M\phi$ is considered, it must first be transformed into the standard form before the HQRI solution technique can be used. As pointed out in Section 10.2.5, this transformation is effective in only some cases.
|
|
|
|
In the following discussion, we consider the problem $K\phi = \lambda\phi$ , in which K may have zero (and could also have negative) eigenvalues. Therefore, it is not necessary to impose a shift prior to applying the HQRI algorithm in order to solve for only positive eigenvalues (see Section 10.2.3). The name “HQRI solution method” stands for the following three solution steps:
|
|
|
|
1. Householder transformations are employed to reduce the matrix K to tridiagonal form.
|
|
2. QR iteration yields all eigenvalues.
|
|
3. Using inverse iteration the required eigenvectors of the tridiagonal matrix are calculated. These vectors are transformed to obtain the eigenvectors of K.
|
|
|
|
A basic difference from the Jacobi solution method is that the matrix is first transformed without iteration into a tridiagonal form. This matrix can then be used effectively in the QR iterative solution, in which all eigenvalues are calculated. Finally, only those eigenvectors that are actually requested are evaluated. We will note that unless many eigenvectors must be calculated, the transformation of K into tridiagonal form requires most of the numerical operations. In the following we consider in detail the three distinct steps carried out in the HQRI solution.
|
|
|
|
The Householder reduction. The Householder reduction to tridiagonal form involves $n - 2$ transformations of the form (11.72); i.e., using $\mathbf{K}_1 = \mathbf{K}$ , we calculate
|
|
|
|
$$
|
|
\mathbf {K} _ {k + 1} = \mathbf {P} _ {k} ^ {T} \mathbf {K} _ {k} \mathbf {P} _ {k}; \quad k = 1, \dots , n - 2 \tag {11.95}
|
|
$$
|
|
|
|
<!-- source-page: 945 -->
|
|
|
|
where the $P_{k}$ are Householder transformation matrices (reflection matrices, see Exercise 2.6):
|
|
|
|
$$
|
|
\mathbf {P} _ {k} = \mathbf {I} - \theta \mathbf {w} _ {k} \mathbf {w} _ {k} ^ {T} \tag {11.96}
|
|
$$
|
|
|
|
$$
|
|
\theta = \frac {2}{\mathbf {w} _ {k} ^ {T} \mathbf {w} _ {k}} \tag {11.97}
|
|
$$
|
|
|
|
To show how the vector $w_{k}$ that defines the matrix $P_{k}$ is calculated, we consider the case k = 1, which is typical. First, we partition $K_{1}, P_{1}$ , and $w_{1}$ into submatrices as follows:
|
|
|
|
$$
|
|
\mathbf {P} _ {1} = \left[ \begin{array}{l l} 1 & 0 \\ 0 & \overline {{\mathbf {P}}} _ {1} \end{array} \right]; \quad \mathbf {w} _ {1} = \left[ \begin{array}{l} 0 \\ \overline {{\mathbf {w}}} _ {1} \end{array} \right] \tag {11.98}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {K} _ {1} = \left[ \begin{array}{l l} k _ {1 1} & \mathbf {k} _ {1} ^ {T} \\ \mathbf {k} _ {1} & \mathbf {K} _ {1 1} \end{array} \right]
|
|
$$
|
|
|
|
where $\mathbf{K}_{11},\overline{\mathbf{P}}_1$ , and $\overline{\mathbf{w}}_1$ are of order $n - 1$ . In the general case of step $k$ , we have corresponding matrices of order $n - k$ . Performing the multiplications in (11.95), we obtain, using the notation of (11.98),
|
|
|
|
$$
|
|
\mathbf {K} _ {2} = \left[ \begin{array}{c c} k _ {1 1} & \mathbf {k} _ {1} ^ {T} \overline {{\mathbf {P}}} _ {1} \\ \overline {{\mathbf {P}}} _ {1} ^ {T} \mathbf {k} _ {1} & \overline {{\mathbf {P}}} _ {1} ^ {T} \mathbf {K} _ {1 1} \overline {{\mathbf {P}}} _ {1} \end{array} \right] \tag {11.99}
|
|
$$
|
|
|
|
The condition is now that the first column and row of $\mathbf{K}_2$ should be in tridiagonal form; i.e., we want $\mathbf{K}_2$ to be in the form
|
|
|
|
$$
|
|
\mathbf {K} _ {2} = \left[ \begin{array}{c c c c c c} k _ {1 1} & \times & 0 & \dots & 0 \\ \hline \times & & & & \\ 0 & & & & \\ \cdot & & & & \\ \cdot & & & \overline {{\mathbf {K}}} _ {2} & \\ 0 & & & & \end{array} \right] \tag {11.100}
|
|
$$
|
|
|
|
where $\times$ indicates a nonzero value and
|
|
|
|
$$
|
|
\overline {{{\mathbf {K}}}} _ {2} = \overline {{{\mathbf {P}}}} _ {1} ^ {T} \mathbf {K} _ {1 1} \overline {{{\mathbf {P}}}} _ {1} \tag {11.101}
|
|
$$
|
|
|
|
The form of $K_{2}$ in (11.100) is achieved by realizing that $\overline{P}_{1}$ is a reflection matrix. Therefore, we can use $\overline{P}_{1}$ to reflect the vector $k_{1}$ of $K_{1}$ in (11.98) into a vector that has only its first component nonzero. Since the length of the new vector must be the length of $k_{1}$ , we determine $\overline{w}_{1}$ from the condition
|
|
|
|
$$
|
|
(\mathbf {I} - \theta \overline {{\mathbf {w}}} _ {1} \overline {{\mathbf {w}}} _ {1} ^ {T}) \mathbf {k} _ {1} = \pm \| \mathbf {k} _ {1} \| _ {2} \mathbf {e} _ {1} \tag {11.102}
|
|
$$
|
|
|
|
where $e_{1}$ is a unit vector of dimension n - 1; i.e., $e_{1}^{T} = [1 \quad 0 \quad 0 \quad \ldots \quad 0]$ , and the + or - sign can be selected to obtain the best numerical stability. Noting that we only need to solve for a multiple of $\overline{w}_{1}$ (i.e., only the direction of the vector normal to the plane of reflection is important, see Exercise 2.6), we obtain from (11.102) as a suitable value for $\overline{w}_{1}$ ,
|
|
|
|
$$
|
|
\overline {{{\mathbf {w}}}} _ {1} = \mathbf {k} _ {1} + \operatorname{sign} \left(k _ {2 1}\right) \| \mathbf {k} _ {1} \| _ {2} \mathbf {e} _ {1} \tag {11.103}
|
|
$$
|
|
|
|
where $k_{21}$ is element (2, 1) of $\mathbf{K}_1$ .
|
|
|
|
With $\overline{w}_{1}$ defined in (11.103), the first Householder transformation, k = 1 in (11.95), can be carried out. In the next step, k = 2, we can consider the matrix $\overline{K}_{2}$ in (11.100) in the same way as we considered $K_{1}$ in (11.98) to (11.103) because the reduction of the first
|
|
|
|
<!-- source-page: 946 -->
|
|
|
|
column and row of $\overline{K}_{2}$ does not affect the first column and row in $K_{2}$ . Thus, the general algorithm for the transformation of K into tridiagonal form is established. We demonstrate the procedure in the following example.
|
|
|
|
EXAMPLE 11.13: Use Householder transformation matrices to reduce K to tridiagonal form, where
|
|
|
|
$$
|
|
\mathbf {K} = \left[ \begin{array}{r r r r} 5 & - 4 & 1 & 0 \\ - 4 & 6 & - 4 & 1 \\ 1 & - 4 & 6 & - 4 \\ 0 & 1 & - 4 & 5 \end{array} \right]
|
|
$$
|
|
|
|
Here, using (11.95) to (11.103) to reduce column 1, we have
|
|
|
|
$$
|
|
\overline {{{\mathbf {w}}}} _ {1} = \left[ \begin{array}{c} - 4 \\ 1 \\ 0 \end{array} \right] - 4. 1 2 3 1 \left[ \begin{array}{l} 1 \\ 0 \\ 0 \end{array} \right] = \left[ \begin{array}{c} - 8. 1 2 3 1 \\ 1 \\ 0 \end{array} \right]
|
|
$$
|
|
|
|
Hence, $\mathbf{w}_1 = \begin{bmatrix} 0 \\ -8.1231 \\ 1 \\ 0 \end{bmatrix}$ ; $\theta_1 = 0.0298575$
|
|
|
|
$$
|
|
\mathbf {P} _ {1} = \left[ \begin{array}{c c c c} 1 & 0 & 0 & 0 \\ 0 & - 0. 9 7 0 1 & 0. 2 4 2 5 & 0 \\ 0 & 0. 2 4 2 5 & 0. 9 7 0 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right]
|
|
$$
|
|
|
|
and $\mathbf{K}_2 = \begin{bmatrix} 5 & 4.1231 & 0 & 0 \\ 4.1231 & 7.8823 & 3.5294 & -1.9403 \\ 0 & 3.5294 & 4.1177 & -3.6380 \\ 0 & -1.9403 & -3.6380 & 5 \end{bmatrix}$
|
|
|
|
Next we reduce column 2,
|
|
|
|
$$
|
|
\overline {{{\mathbf {w}}}} _ {2} = \left[ \begin{array}{c} 3. 5 2 9 4 \\ - 1. 9 4 0 3 \end{array} \right] + 4. 0 2 7 6 \left[ \begin{array}{l} 1 \\ 0 \end{array} \right] = \left[ \begin{array}{c} 7. 5 5 7 0 \\ - 1. 9 4 0 3 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\mathbf {w} _ {2} = \left[ \begin{array}{c} 0 \\ 0 \\ 7. 5 5 7 0 \\ - 1. 9 4 0 3 \end{array} \right]; \quad \theta_ {2} = 0. 0 3 2 8 5 5 3
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {2} = \left[ \begin{array}{c c c c} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & - 0. 8 7 6 3 & 0. 4 8 1 7 \\ 0 & 0 & 0. 4 8 1 7 & 0. 8 7 6 3 \end{array} \right]
|
|
$$
|
|
|
|
Hence, $\mathbf{K}_3 = \begin{bmatrix} 5 & 4.1231 & 0 & 0 \\ 4.1231 & 7.8823 & -4.0276 & 0 \\ 0 & -4.0276 & 7.3941 & 2.3219 \\ 0 & 0 & 2.3219 & 1.7236 \end{bmatrix}$
|
|
|
|
<!-- source-page: 947 -->
|
|
|
|
Some important numerical aspects should be noted. First, the reduced matrices $K_{2}$ , $K_{3}$ , $\ldots$ , $K_{n-1}$ are symmetric. This means that in the reduction we need to store only the lower symmetric part of K. Furthermore, to store the $\overline{w}_{k}$ , $k = 1, 2, \ldots, n - 2$ , we can use the storage locations below the subdiagonal elements in the matrix currently being reduced.
|
|
|
|
A disadvantage of the Householder transformations is that the bandwidth is increased in the unreduced part of $K_{k+1}$ . Hence, in the reduction, essentially no advantage can be taken of the bandedness of K.
|
|
|
|
An important aspect of the transformation is the evaluation of the matrix product $\overline{P}_{1}^{T}K_{11}\overline{P}_{1}$ and the similar products required in the next steps. In the most general case a triple matrix product with matrices of order n requires $2n^{3}$ operations, and if this many operations were required, the Householder reduction would be quite uneconomical. However, by taking advantage of the special nature of the matrix $\overline{P}_{1}$ , we can evaluate the product $\overline{P}_{1}^{T}K_{11}\overline{P}_{1}$ by calculating
|
|
|
|
$$
|
|
\left. \begin{array}{l} \mathbf {v} _ {1} = \mathbf {K} _ {1 1} \overline {{\mathbf {w}}} _ {1} \\ \mathbf {p} _ {1} ^ {T} = \theta_ {1} \mathbf {v} _ {1} ^ {T} \\ \beta_ {1} = \mathbf {p} _ {1} ^ {T} \overline {{\mathbf {w}}} _ {1} \\ \mathbf {q} _ {1} = \mathbf {p} _ {1} - \theta_ {1} \beta_ {1} \overline {{\mathbf {w}}} _ {1} \end{array} \right\} \tag {11.104}
|
|
$$
|
|
|
|
and then $\overline{\mathbf{P}}_1^T\mathbf{K}_{11}\overline{\mathbf{P}}_1 = \mathbf{K}_{11} - \overline{\mathbf{w}}_1\mathbf{p}_1^T -\mathbf{q}_1\overline{\mathbf{w}}_1^T$ (11.105)
|
|
|
|
which requires only about $3m^{2} + 3m$ operations, where $m$ is the order of $\overline{\mathbf{P}}_{1}$ and $\mathbf{K}_{11}$ (i.e., $m = n - 1$ in this case). Hence, the multiplication $\overline{\mathbf{P}}_{1}^{7}[\mathbf{K}_{11}\overline{\mathbf{P}}_{1}$ requires a number of operations of the order $m$ -squared rather than $m$ -cubed, which is a very significant reduction. We demonstrate the procedure given in (11.104) and (11.105) by reworking Example 11.13.
|
|
|
|
EXAMPLE 11.14: Use the relations given in (11.104) and (11.105) to reduce the matrix K in Example 11.13 to tridiagonal form.
|
|
|
|
Here we obtain for the reduction of column 1 using $\overline{w}_{1}$ and $\theta_{1}$ calculated in Example 11.13,
|
|
|
|
$$
|
|
\begin{array}{l} \mathbf {v} _ {1} = \left[ \begin{array}{r r r} 6 & - 4 & 1 \\ - 4 & 6 & - 4 \\ 1 & - 4 & 5 \end{array} \right] \overline {{\mathbf {w}}} _ {1} = \left[ \begin{array}{r} - 5 2. 7 3 8 \\ 3 8. 4 9 2 4 \\ - 1 2. 1 2 3 1 \end{array} \right] \\ \mathbf {p} _ {1} ^ {T} = \left[ \begin{array}{l l l} - 1. 5 7 4 6 & 1. 1 4 9 3 & - 0. 3 6 1 9 7 \end{array} \right]; \quad \beta_ {1} = 1 3. 9 4 0 3; \quad \mathbf {q} _ {1} = \left[ \begin{array}{c} 1. 8 0 6 4 \\ 0. 7 3 3 1 \\ - 0. 3 6 2 0 \end{array} \right] \\ \end{array}
|
|
$$
|
|
|
|
$$
|
|
\begin{array}{l} \overline {{{\mathbf {P}}}} _ {1} ^ {T} \mathbf {K} _ {1 1} \overline {{{\mathbf {P}}}} _ {1} = \left[ \begin{array}{c c c} 6 & - 4 & 1 \\ - 4 & 6 & - 4 \\ 1 & - 4 & 5 \end{array} \right] - \left[ \begin{array}{c c c} 1 2. 7 9 1 0 & - 9. 3 3 5 & 2. 9 4 0 3 \\ - 1. 5 7 4 6 & 1. 1 4 9 3 & - 0. 3 6 2 0 \\ 0 & 0 & 0 \end{array} \right] \\ - \left[ \begin{array}{c c c} - 1 4. 6 7 3 4 & 1. 8 0 6 4 & 0 \\ - 5. 9 5 4 8 & 0. 7 3 3 0 7 & 0 \\ 2. 9 4 0 3 & - 0. 3 6 2 0 & 0 \end{array} \right] \\ \end{array}
|
|
$$
|
|
|
|
$$
|
|
\overline {{{\mathbf {P}}}} _ {1} ^ {T} \mathbf {K} _ {1 1} \overline {{{\mathbf {P}}}} _ {1} = \left[ \begin{array}{c c c} 7. 8 8 2 3 & 3. 5 2 9 4 & - 1. 9 4 0 3 \\ 3. 5 2 9 4 & 4. 1 1 7 7 & - 3. 6 3 8 0 \\ - 1. 9 4 0 3 & - 3. 6 3 8 0 & 5 \end{array} \right]
|
|
$$
|
|
|
|
<!-- source-page: 948 -->
|
|
|
|
and hence, $\mathbf{K}_2 = \begin{bmatrix} 5 & 4.1231 & 0 & 0 \\ 4.1231 & 7.8823 & 3.5294 & -1.9403 \\ 0 & 3.5294 & 4.1177 & -3.6380 \\ 0 & -1.9403 & -3.6380 & 5 \end{bmatrix}$
|
|
|
|
Next we reduce the second column
|
|
|
|
$$
|
|
\begin{array}{l} \mathbf {v} _ {2} = \left[ \begin{array}{c c} 4. 1 1 7 7 & - 3. 6 3 8 0 \\ - 3. 6 3 8 0 & 5 \end{array} \right] \overline {{\mathbf {w}}} _ {2} = \left[ \begin{array}{c} 3 8. 1 7 5 9 \\ - 3 7. 1 9 4 1 \end{array} \right] \\ \mathbf {p} _ {2} ^ {T} = [ 1. 2 5 4 3 - 1. 2 2 2 0 ]; \quad \beta_ {2} = 1 1. 8 4 9 7; \quad \mathbf {q} _ {2} = \left[ \begin{array}{l} - 1. 6 8 7 8 \\ - 0. 4 6 6 6 \end{array} \right] \\ \end{array}
|
|
$$
|
|
|
|
$$
|
|
\begin{array}{l} \overline {{{\mathbf {P}}}} _ {2} ^ {T} \mathbf {K} _ {2 2} \overline {{{\mathbf {P}}}} _ {2} = \left[ \begin{array}{c c} 4. 1 1 7 7 & - 3. 6 3 8 0 \\ - 3. 6 3 8 0 & 5 \end{array} \right] - \left[ \begin{array}{c c} 9. 4 7 8 6 & - 9. 2 3 4 8 \\ - 2. 4 3 3 7 & 2. 3 7 1 1 \end{array} \right] - \left[ \begin{array}{c c} - 1 2. 7 5 5 0 & 3. 2 7 4 9 \\ - 3. 5 2 6 3 & 0. 9 0 5 3 8 \end{array} \right] \\ = \left[ \begin{array}{l l} 7. 3 9 4 1 & 2. 3 2 1 9 \\ 2. 3 2 1 9 & 1. 7 2 3 6 \end{array} \right] \\ \end{array}
|
|
$$
|
|
|
|
and hence, $\mathbf{K}_3 = \begin{bmatrix} 5 & 4.1231 & 0 & 0 \\ 4.1231 & 7.9823 & -4.0276 & 0 \\ 0 & -4.0276 & 7.3941 & 2.3219 \\ 0 & 0 & 2.3219 & 1.7236 \end{bmatrix}$
|
|
|
|
The QR iteration. In the HQRI solution procedure, the QR iteration is applied to the tridiagonal matrix obtained by the Householder transformation of K. However, it should be realized that the QR iteration could be applied to the original matrix K as well, and that the transformation of K into tridiagonal form prior to the iteration is merely carried out to improve the efficiency of solution. In the following we therefore consider first how the iteration is applied to a general symmetric matrix K.
|
|
|
|
The name "QR iteration" derives from the notation used in the algorithm. Namely, the basic step in the iteration is to decompose $\mathbf{K}$ into the form
|
|
|
|
$$
|
|
\mathbf {K} = \mathbf {Q R} \tag {11.106}
|
|
$$
|
|
|
|
where $\mathbf{Q}$ is an orthogonal and $\mathbf{R}$ is an upper triangular matrix. We then form
|
|
|
|
$$
|
|
\mathbf {R Q} = \mathbf {Q} ^ {T} \mathbf {K Q} \tag {11.107}
|
|
$$
|
|
|
|
Therefore, by calculating RQ we in fact carry out a transformation of the form (11.72).
|
|
|
|
The factorization in (11.106) could be obtained by applying the Gram-Schmidt process to the columns of K. In practice, it is more effective to reduce K into upper triangular form using Jacobi rotation matrices; i.e., we evaluate
|
|
|
|
$$
|
|
\mathbf {P} _ {n, n - 1} ^ {T} \dots \mathbf {P} _ {3, 1} ^ {T} \mathbf {P} _ {2, 1} ^ {T} \mathbf {K} = \mathbf {R} \tag {11.108}
|
|
$$
|
|
|
|
where the rotation matrix $\mathbf{P}_{j,i}^{T}$ is selected to zero element $(j,i)$ . Using (11.108), we have, corresponding to (11.106),
|
|
|
|
$$
|
|
\mathbf {Q} = \mathbf {P} _ {2, 1} \mathbf {P} _ {3, 1} \dots \mathbf {P} _ {n, n - 1} \tag {11.109}
|
|
$$
|
|
|
|
The QR iteration algorithm is obtained by repeating the process given in (11.106) and (11.107). Using the notation $\mathbf{K}_1 = \mathbf{K}$ , we form
|
|
|
|
$$
|
|
\mathbf {K} _ {k} = \mathbf {Q} _ {k} \mathbf {R} _ {k} \tag {11.110}
|
|
$$
|
|
|
|
<!-- source-page: 949 -->
|
|
|
|
and then
|
|
|
|
$$
|
|
\mathbf {K} _ {k + 1} = \mathbf {R} _ {k} \mathbf {Q} _ {k} \tag {11.111}
|
|
$$
|
|
|
|
where, then, disregarding that eigenvalues and eigenvectors may not be in the usual order,
|
|
|
|
$$
|
|
\mathbf {K} _ {k + 1} \rightarrow \Lambda \quad \text { and } \quad \mathbf {Q} _ {1} \dots \mathbf {Q} _ {k - 1} \mathbf {Q} _ {k} \rightarrow \Phi \quad \text { as } k \rightarrow \infty
|
|
$$
|
|
|
|
We demonstrate the iteration process in the following example.
|
|
|
|
EXAMPLE 11.15: Use the QR iteration with Q obtained as a product of Jacobi rotation matrices to calculate the eigensystem of K, where
|
|
|
|
$$
|
|
\mathbf {K} = \left[ \begin{array}{r r r r} 5 & - 4 & 1 & 0 \\ - 4 & 6 & - 4 & 1 \\ 1 & - 4 & 6 & - 4 \\ 0 & 1 & - 4 & 5 \end{array} \right]
|
|
$$
|
|
|
|
The Jacobi rotation matrix $P_{j,i}^{T}$ to zero element $(j, i)$ in the current matrix is given by
|
|
|
|
$$
|
|
\mathbf {P} _ {j, i} ^ {T} = \left[ \begin{array}{c c c c c c} 1 & & & & & \\ & \ddots & & & & \\ & & 1 & & & \\ & & \cos \theta & \sin \theta & - & - i \text {th row} \\ & & & 1 & & \\ & & & \ddots & & \\ & & & & 1 & \\ & - \sin \theta & & \cos \theta & - & - j \text {th row} \\ & & & & 1 & \\ & & & & & \ddots \\ & & & & & 1 \end{array} \right]
|
|
$$
|
|
|
|
where $\sin \theta = \frac{\overline{k}_{ji}}{(\overline{k}_{ii}^2 + \overline{k}_{ji}^2)^{1 / 2}};$ $\cos \theta = \frac{\overline{k}_{ii}}{(\overline{k}_{ii}^2 + \overline{k}_{ji}^2)^{1 / 2}}$
|
|
|
|
and the bar indicates that the elements of the current matrix are used.
|
|
|
|
Proceeding as in (11.108), we obtain for element (2, 1),
|
|
|
|
$$
|
|
\sin \theta = - 0. 6 2 4 7; \quad \cos \theta = 0. 7 8 0 9
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {2, 1} = \left[ \begin{array}{c c c c} 0. 7 8 0 9 & 0. 6 2 4 7 & 0 & 0 \\ - 0. 6 2 4 7 & 0. 7 8 0 9 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right]
|
|
$$
|
|
|
|
and $\mathbf{P}_{2,1}^{T}\mathbf{K} = \begin{bmatrix} 6.403 & -6.872 & 3.280 & -0.6247\\ 0 & 2.186 & -2.499 & 0.7809\\ 1 & -4 & 6 & -4\\ 0 & 1 & -4 & 5 \end{bmatrix}$
|
|
|
|
Next we zero element (3, 1),
|
|
|
|
$$
|
|
\mathbf {P} _ {3, 1} ^ {T} \mathbf {P} _ {2, 1} ^ {T} \mathbf {K} = \left[ \begin{array}{c c c c} 6. 4 8 1 & - 7. 4 0 7 & 4. 1 6 6 & - 1. 2 3 4 \\ 0 & 2. 1 8 6 & - 2. 4 9 9 & 0. 7 8 0 9 \\ 0 & - 2. 8 9 2 & 5. 4 2 2 & - 3. 8 5 6 \\ 0 & 1 & - 4 & 5 \end{array} \right]
|
|
$$
|
|
|
|
<!-- source-page: 950 -->
|
|
|
|
Noting that element (4, 1) is zero already, we continue with the factorization by zeroing element (3, 2),
|
|
|
|
$$
|
|
\mathbf {P} _ {3, 2} ^ {T} \mathbf {P} _ {3, 1} ^ {T} \mathbf {P} _ {2, 1} ^ {T} \mathbf {K} = \left[ \begin{array}{c c c c} 6. 4 8 1 & - 7. 4 0 7 & 4. 1 6 6 & - 1. 2 3 4 \\ 0 & 3. 6 2 5 & - 5. 8 3 2 & 3. 5 4 6 \\ 0 & 0 & 1. 2 7 7 & - 1. 7 0 3 \\ 0 & 1 & - 4 & 5 \end{array} \right]
|
|
$$
|
|
|
|
Proceeding in the same manner, we obtain
|
|
|
|
$$
|
|
\mathbf {P} _ {4, 2} ^ {T} \mathbf {P} _ {3, 2} ^ {T} \mathbf {P} _ {3, 1} ^ {T} \mathbf {P} _ {2, 1} ^ {T} \mathbf {K} = \left[ \begin{array}{c c c c} 6. 4 8 1 & - 7. 4 0 7 & 4. 1 6 6 & - 1. 2 3 4 \\ 0 & 3. 7 6 1 & - 6. 6 8 6 & 4. 7 4 8 \\ 0 & 0 & 1. 2 7 7 & - 1. 7 0 3 \\ 0 & 0 & - 2. 3 0 5 & 3. 8 7 7 \end{array} \right]
|
|
$$
|
|
|
|
and finally, $\mathbf{R}_1 = \mathbf{P}_{4,3}^T\mathbf{P}_{4,2}^T\mathbf{P}_{3,2}^T\mathbf{P}_{3,1}^T\mathbf{P}_{2,1}^T\mathbf{K}$
|
|
|
|
$$
|
|
\mathbf {R} _ {1} = \left[ \begin{array}{c c c c} 6. 4 8 1 & - 7. 4 0 7 & 4. 1 6 6 & - 1. 2 3 4 \\ 0 & 3. 7 6 1 & - 6. 6 8 6 & 4. 7 4 8 \\ 0 & 0 & 2. 6 3 5 & - 4. 2 1 6 \\ 0 & 0 & 0 & 0. 3 8 9 2 \end{array} \right]
|
|
$$
|
|
|
|
Also, we have $\mathbf{Q}_1 = \mathbf{P}_{2,1}\mathbf{P}_{3,1}\mathbf{P}_{3,2}\mathbf{P}_{4,2}\mathbf{P}_{4,3}$
|
|
|
|
$$
|
|
\mathbf {Q} _ {1} = \left[ \begin{array}{c c c c} 0. 7 7 1 5 & 0. 4 5 5 8 & 0. 3 1 6 2 & 0. 3 1 1 4 \\ - 0. 6 1 7 2 & 0. 3 7 9 9 & 0. 4 2 1 6 & 0. 5 4 4 9 \\ 0. 1 5 4 3 & - 0. 7 5 9 7 & 0. 1 0 5 4 & 0. 6 2 2 8 \\ 0 & 0. 2 6 5 9 & - 0. 8 4 3 3 & 0. 4 6 7 1 \end{array} \right]
|
|
$$
|
|
|
|
The first iteration step of QR is completed by calculating
|
|
|
|
$$
|
|
\mathbf {K} _ {2} = \mathbf {R} _ {1} \mathbf {Q} _ {1}
|
|
$$
|
|
|
|
to obtain $\mathbf{K}_2 = \begin{bmatrix} 10.21 & -3.353 & 0.4066 & 0 \\ -3.353 & 7.771 & -3.123 & 0.1035 \\ 0.4066 & -3.123 & 3.833 & -0.3282 \\ 0 & 0.1035 & -0.3282 & 0.1818 \end{bmatrix}$
|
|
|
|
The following results are obtained in the next steps of QR. For k = 2:
|
|
|
|
$$
|
|
\mathbf {R} _ {2} = \left[ \begin{array}{c c c c} 1 0. 7 6 & - 5. 7 2 3 & 1. 5 0 4 & - 0. 0 4 4 6 \\ 0 & 6. 9 7 4 & - 4. 1 6 3 & 0. 2 2 8 4 \\ 0 & 0 & 2. 2 6 5 & - 0. 2 7 5 2 \\ 0 & 0 & 0 & 0. 1 4 7 1 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\mathbf {Q} _ {1} \mathbf {Q} _ {2} = \left[ \begin{array}{c c c c} 0. 6 0 2 4 & 0. 4 9 4 3 & 0. 5 0 8 4 & 0. 3 6 6 5 \\ - 0. 6 8 8 5 & - 0. 0 2 5 7 & 0. 4 0 9 9 & 0. 5 9 7 8 \\ 0. 3 8 7 3 & - 0. 6 4 0 9 & - 0. 2 7 1 5 & 0. 6 0 4 6 \\ - 0. 1 1 4 7 & 0. 5 8 6 7 & 0. 7 0 7 0 & 0. 3 7 7 9 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\mathbf {K} _ {3} = \mathbf {R} _ {2} \mathbf {Q} _ {2} = \left[ \begin{array}{c c c c} 1 2. 0 5 & - 2. 3 3 1 & 0. 0 8 5 6 & 0 \\ - 2. 3 3 1 & 7. 7 2 6 & - 0. 9 4 8 3 & 0. 0 0 2 2 \\ 0. 0 8 5 6 & - 0. 9 4 8 3 & 2. 0 7 4 0 & - 0. 0 1 7 3 \\ 0 & 0. 0 0 2 2 & - 0. 0 1 7 3 & 0. 1 4 6 1 \end{array} \right]
|
|
$$
|