김경종
421 lines
28 KiB
Markdown
421 lines
28 KiB
Markdown
<!-- source-page: 931 -->
|
|
|
|
A procedure that has been used very effectively is the threshold Jacobi method, in which the off-diagonal elements are tested sequentially, namely, row by row (or column by column), and a rotation is applied only if the element is larger than the threshold for that sweep. To define an appropriate threshold we note that, physically, in the diagonalization of K we want to reduce the coupling between the degrees of freedom i and j. A measure of this coupling is given by $(k_{ij}^{2}/k_{ii}k_{jj})^{1/2}$ , and it is this factor that can be used effectively in deciding whether to apply a rotation. In addition to having a realistic threshold tolerance, it is also necessary to measure convergence. As described above, $K_{k+1} \to \Lambda$ as $k \to \infty$ , but in the numerical computations we seek only a close enough approximation to the eigenvalues and corresponding eigenvectors. Let l be the last iteration; i.e., we have, to the precision required,
|
|
|
|
$$
|
|
\mathbf {K} _ {l + 1} \doteq \mathbf {\Lambda} \tag {11.82}
|
|
$$
|
|
|
|
Then we say that convergence to a tolerance $s$ has been achieved provided that
|
|
|
|
$$
|
|
\frac {\left| k _ {i i} ^ {(l + 1)} - k _ {i i} ^ {(l)} \right|}{k _ {i i} ^ {(l + 1)}} \leq 1 0 ^ {- s}; \quad i = 1, \dots , n \tag {11.83}
|
|
$$
|
|
|
|
and
|
|
|
|
$$
|
|
\left[ \frac {(k _ {i j} ^ {(l + 1)}) ^ {2}}{k _ {i i} ^ {(l + 1)} k _ {j j} ^ {(l + 1)}} \right] ^ {1 / 2} \leq 1 0 ^ {- s}; \quad \text { all } i, j; i < j \tag {11.84}
|
|
$$
|
|
|
|
The relation in (11.83) has to be satisfied because the element $k_{ii}^{(l+1)}$ is the current approximation to an eigenvalue, and the relation states that the current and last approximations to the eigenvalues did not change in the first s digits. This convergence measure is essentially the same as the one used in vector iteration in (11.20). The relation in (11.84) ensures that the off-diagonal elements are indeed small.
|
|
|
|
Having discussed the main aspects of the iteration, we may now summarize the actual solution procedure. The following steps have been used in a threshold Jacobi iteration.
|
|
|
|
1. Initialize the threshold for the sweep. Typically, the threshold used for sweep $m$ may be $10^{-2m}$ .
|
|
2. For all $i, j$ with $i < j$ calculate the coupling factor $[(k_{ij}^{(k)})^2 / k_{ii}^{(k)} k_{jj}^{(k)}]^{1/2}$ and apply a transformation if the factor is larger than the current threshold.
|
|
3. Use (11.83) to check for convergence. If the relation in (11.83) is not satisfied, continue with the next sweep; i.e., go to step 1. If (11.83) is satisfied, check if (11.84) is also satisfied; if "yes," the iteration converged; if "no," continue with the next sweep.
|
|
|
|
So far we have stated the algorithm but we have not shown that convergence will indeed always occur. The proof of convergence has been given elsewhere (see J. H. Wilkinson [A]) and will not be repeated here because little additional insight into the working of the solution procedure would be gained. However, one important point should be noted—that convergence is quadratic once the off-diagonal elements are small. Since rapid convergence is obtained once the off-diagonal elements are small, little extra cost is involved in solving for the eigensystem to high accuracy when an approximate solution has been obtained. In practical solutions we use m = 2 and s = 12, and about six sweeps are required for solution of the eigensystem to high accuracy. A program used is given in the next section, when we discuss the solution of the generalized eigenproblem $K\phi = \lambda M\phi$ .
|
|
|
|
<!-- source-page: 932 -->
|
|
|
|
EXAMPLE 11.9: Calculate the eigensystem of the matrix 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]
|
|
$$
|
|
|
|
Use the threshold Jacobi iteration described above.
|
|
|
|
To demonstrate the solution algorithm we give one sweep in detail and then the results obtained in the next sweeps.
|
|
|
|
For sweep 1 we have as a threshold $10^{-2}$ . We therefore obtain the following results. For $i = 1, j = 2$ :
|
|
|
|
$$
|
|
\cos \theta = 0. 7 4 9 7; \quad \sin \theta = 0. 6 6 1 8
|
|
$$
|
|
|
|
and thus
|
|
|
|
$$
|
|
\mathbf {P} _ {1} = \left[ \begin{array}{c c c c} 0. 7 4 9 7 & - 0. 6 6 1 8 & 0 & 0 \\ 0. 6 6 1 8 & 0. 7 4 9 7 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {1} ^ {\dagger} \mathbf {K P} _ {1} = \left[ \begin{array}{c c c c} 1. 4 6 9 & 0 & - 1. 8 9 8 & 0. 6 6 1 8 \\ 0 & 9. 5 3 1 & - 3. 6 6 1 & 0. 7 4 9 7 \\ - 1. 8 9 8 & - 3. 6 6 1 & 6 & - 4 \\ 0. 6 6 1 8 & 0. 7 4 9 7 & - 4 & 5 \end{array} \right]
|
|
$$
|
|
|
|
For $i = 1, j = 3$ :
|
|
|
|
$$
|
|
\cos \theta = 0. 9 3 9 8; \quad \sin \theta = 0. 3 4 1 6
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {2} = \left[ \begin{array}{c c c c} 0. 9 3 9 8 & 0 & - 0. 3 4 1 6 & 0 \\ 0 & 1 & 0 & 0 \\ 0. 3 4 1 6 & 0 & 0. 9 3 9 8 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {2} ^ {T} \mathbf {P} _ {1} ^ {T} \mathbf {K P} _ {1} \mathbf {P} _ {2} = \left[ \begin{array}{c c c c} 0. 7 7 9 2 & - 1. 2 5 0 & 0 & - 0. 7 4 4 4 \\ - 1. 2 5 0 & 9. 5 3 1 & - 3. 4 4 0 & 0. 7 4 9 7 \\ 0 & - 3. 4 4 0 & 6. 6 9 0 & - 3. 9 8 6 \\ - 0. 7 4 4 4 & 0. 7 4 9 7 & - 3. 9 8 6 & 5 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {1} \mathbf {P} _ {2} = \left[ \begin{array}{c c c c} 0. 7 0 4 6 & - 0. 6 6 1 8 & - 0. 2 5 6 1 & 0 \\ 0. 6 2 2 0 & 0. 7 4 9 7 & - 0. 2 2 6 1 & 0 \\ 0. 3 4 1 6 & 0 & 0. 9 3 9 8 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right]
|
|
$$
|
|
|
|
For $i = 1, j = 4$ :
|
|
|
|
$$
|
|
\cos \theta = 0. 9 8 5 7; \quad \sin \theta = 0. 1 6 8 7
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {3} = \left[ \begin{array}{c c c c} 0. 9 8 5 7 & 0 & 0 & - 0. 1 6 8 7 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0. 1 6 8 7 & 0 & 0 & 0. 9 8 5 7 \end{array} \right]
|
|
$$
|
|
|
|
<!-- source-page: 933 -->
|
|
|
|
$$
|
|
\mathbf {P} _ {3} ^ {T} \mathbf {P} _ {2} ^ {T} \mathbf {P} _ {1} ^ {T} \mathbf {K P} _ {1} \mathbf {P} _ {2} \mathbf {P} _ {3} = \left[ \begin{array}{c c c c} 0. 6 5 1 8 & - 1. 1 0 6 & - 0. 6 7 2 5 & 0 \\ - 1. 1 0 6 & 9. 5 3 1 & - 3. 4 4 0 & 0. 9 4 9 9 \\ - 0. 6 7 2 5 & - 3. 4 4 0 & 6. 6 9 0 & - 3. 9 2 8 \\ 0 & 0. 9 4 9 9 & - 3. 9 2 8 & 5. 1 2 7 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {1} \mathbf {P} _ {2} \mathbf {P} _ {3} = \left[ \begin{array}{c c c c} 0. 6 9 4 5 & - 0. 6 6 1 8 & - 0. 2 5 6 1 & - 0. 1 1 8 9 \\ 0. 6 1 3 1 & 0. 7 4 9 7 & - 0. 2 2 6 1 & - 0. 1 0 5 0 \\ 0. 3 3 6 7 & 0 & 0. 9 3 9 8 & - 0. 0 5 7 6 \\ 0. 1 6 8 7 & 0 & 0 & 0. 9 8 5 7 \end{array} \right]
|
|
$$
|
|
|
|
For $i = 2, j = 3$ :
|
|
|
|
$$
|
|
\cos \theta = 0. 8 3 1 2; \quad \sin \theta = - 0. 5 5 6 0
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {4} = \left[ \begin{array}{c c c c} 1 & 0 & 0 & 0 \\ 0 & 0. 8 3 1 2 & 0. 5 5 6 0 & 0 \\ 0 & - 0. 5 5 6 0 & 0. 8 3 1 2 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {4} ^ {T} \mathbf {P} _ {3} ^ {T} \mathbf {P} _ {2} ^ {T} \mathbf {P} _ {1} ^ {T} \mathbf {K P} _ {1} \mathbf {P} _ {2} \mathbf {P} _ {3} \mathbf {P} _ {4} = \left[ \begin{array}{c c c c} 0. 6 5 1 8 & 0. 5 4 5 3 & - 1. 1 7 4 & 0 \\ - 0. 5 4 5 3 & 1 1. 8 3 & 0 & 2. 9 7 4 \\ - 1. 1 7 4 & 0 & 4. 3 8 8 & - 2. 7 3 7 \\ 0 & 2. 9 7 4 & - 2. 7 3 7 & 5. 1 2 7 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {1} \mathbf {P} _ {2} \mathbf {P} _ {3} \mathbf {P} _ {4} = \left[ \begin{array}{c c c c} 0. 6 9 4 5 & - 0. 4 0 7 7 & - 0. 5 8 0 8 & - 0. 1 1 8 9 \\ 0. 6 1 3 1 & 0. 7 4 8 8 & 0. 2 2 8 9 & - 0. 1 0 5 0 \\ 0. 3 3 6 7 & - 0. 5 2 2 6 & 0. 7 8 1 2 & - 0. 0 5 7 6 \\ 0. 1 6 8 2 & 0 & 0 & 0. 9 8 5 7 \end{array} \right]
|
|
$$
|
|
|
|
For $i = 2, j = 4$ :
|
|
|
|
$$
|
|
\cos \theta = 0. 9 3 4 9; \quad \sin \theta = 0. 3 5 4 9
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {5} = \left[ \begin{array}{c c c c} 1 & 0 & 0 & 0 \\ 0 & 0. 9 3 4 9 & 0 & - 0. 3 5 4 9 \\ 0 & 0 & 1 & 0 \\ 0 & 0. 3 5 4 9 & 0 & 0. 9 3 4 9 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {5} ^ {T} \mathbf {P} _ {4} ^ {T} \mathbf {P} _ {3} ^ {T} \mathbf {P} _ {2} ^ {T} \mathbf {P} _ {1} ^ {T} \mathbf {K P} _ {1} \mathbf {P} _ {2} \mathbf {P} _ {3} \mathbf {P} _ {4} \mathbf {P} _ {5} = \left[ \begin{array}{c c c c} 0. 6 5 1 8 & 0. 5 0 9 8 & - 1. 1 7 4 & 0. 1 9 3 5 \\ - 0. 5 0 9 8 & 1 2. 9 6 & 0. 9 7 1 3 & 0 \\ - 1. 1 7 4 & - 0. 9 7 1 3 & 4. 3 8 8 & - 2. 5 5 9 \\ 0. 1 9 3 5 & 0 & - 2. 5 5 9 & 3. 9 9 9 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {1} \mathbf {P} _ {2} \mathbf {P} _ {3} \mathbf {P} _ {4} \mathbf {P} _ {5} = \left[ \begin{array}{c c c c} 0. 6 9 4 5 & - 0. 4 2 3 3 & - 0. 5 8 0 8 & 0. 0 3 3 5 \\ 0. 6 1 3 1 & 0. 6 6 2 8 & 0. 2 2 8 9 & - 0. 3 6 3 9 \\ 0. 3 3 6 7 & 0. 5 0 9 0 & 0. 7 8 1 2 & 0. 1 3 1 6 \\ 0. 1 6 8 7 & 0. 3 4 9 8 & 0 & 0. 9 2 1 3 \end{array} \right]
|
|
$$
|
|
|
|
<!-- source-page: 934 -->
|
|
|
|
To complete the sweep, we zero element (3, 4), using
|
|
|
|
$$
|
|
\cos \theta = 0. 7 3 3 5; \quad \sin \theta = - 0. 6 7 9 7
|
|
$$
|
|
|
|
$$
|
|
\mathbf {P} _ {6} = \left[ \begin{array}{c c c c} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0. 7 3 3 5 & 0. 6 7 9 7 \\ 0 & 0 & - 0. 6 7 9 7 & 0. 7 3 3 5 \end{array} \right]
|
|
$$
|
|
|
|
and hence, the approximations obtained for $\Lambda$ and $\Phi$ are
|
|
|
|
$$
|
|
\mathbf {\Lambda} \doteq \mathbf {P} _ {6} ^ {T} \dots \mathbf {P} _ {1} ^ {T} \mathbf {K} \mathbf {P} _ {1} \dots \mathbf {P} _ {6}
|
|
$$
|
|
|
|
i.e., $\Lambda \doteq \left[ \begin{array}{cccc}0.6518 & -0.5098 & -0.9926 & -0.6560\\ -0.5098 & 12.96 & -0.7124 & -0.6602\\ -0.9926 & -0.7124 & 6.7596 & 0\\ -0.6560 & -0.6602 & 0 & 1.6272 \end{array} \right]$
|
|
|
|
and $\Phi \doteq \mathbf{P}_1\ldots \mathbf{P}_6$
|
|
|
|
i.e., $\Phi \doteq \left[ \begin{array}{cccc}0.6945 & -0.4233 & -0.4488 & -0.3702\\ 0.6131 & 0.6628 & 0.4152 & -0.1113\\ 0.3367 & -0.5090 & 0.4835 & 0.6275\\ 0.1687 & 0.3498 & -0.6264 & 0.6759 \end{array} \right]$
|
|
|
|
After the second sweep we obtain
|
|
|
|
$$
|
|
\boldsymbol {\Lambda} \doteq \left[ \begin{array}{c c c c} 0. 1 5 6 3 & - 0. 3 6 3 5 & 0. 0 0 6 3 & - 0. 0 1 7 6 \\ - 0. 3 6 3 5 & 1 3. 0 8 & - 0. 0 0 2 0 & 0 \\ 0. 0 0 6 3 & - 0. 0 0 2 0 & 6. 8 4 5 & 0 \\ - 0. 0 1 7 6 & 0 & 0 & 1. 9 1 0 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\Phi \doteq \left[ \begin{array}{c c c c} 0. 3 8 7 5 & - 0. 3 6 1 2 & - 0. 6 0 1 7 & - 0. 5 9 7 8 \\ 0. 5 8 8 4 & 0. 6 1 8 4 & 0. 3 7 1 0 & - 0. 3 6 5 7 \\ 0. 6 1 4 8 & - 0. 5 8 4 3 & 0. 3 7 1 4 & 0. 3 7 7 7 \\ 0. 3 5 4 6 & 0. 3 8 1 6 & - 0. 6 0 2 0 & 0. 6 0 5 2 \end{array} \right]
|
|
$$
|
|
|
|
And after the third sweep we have
|
|
|
|
$$
|
|
\boldsymbol {\Lambda} \doteq \left[ \begin{array}{c c c c} 0. 1 4 5 9 & & & \\ & 1 3. 0 9 & & \\ & & 6. 8 5 4 & \\ & & & 1. 9 1 0 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\boldsymbol {\Phi} \doteq \left[ \begin{array}{c c c c} 0. 3 7 1 7 & - 0. 3 7 1 7 & - 0. 6 0 1 5 & - 0. 6 0 1 5 \\ 0. 6 0 1 5 & 0. 6 0 1 5 & 0. 3 7 1 7 & - 0. 3 7 1 7 \\ 0. 6 0 1 5 & - 0. 6 0 1 5 & 0. 3 7 1 7 & 0. 3 7 1 7 \\ 0. 3 7 1 7 & 0. 3 7 1 7 & - 0. 6 0 1 5 & 0. 6 0 1 5 \end{array} \right]
|
|
$$
|
|
|
|
<!-- source-page: 935 -->
|
|
|
|
The approximation for $\Lambda$ is diagonal to the precision given, and we can use
|
|
|
|
$$
|
|
\lambda_ {1} \doteq 0. 1 4 5 9; \quad \boldsymbol {\phi} _ {1} \doteq \left[ \begin{array}{l} 0. 3 7 1 7 \\ 0. 6 0 1 5 \\ 0. 6 0 1 5 \\ 0. 3 7 1 7 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\lambda_ {2} \doteq 1. 9 1 0; \quad \boldsymbol {\Phi} _ {2} \doteq \left[ \begin{array}{c} - 0. 6 0 1 5 \\ - 0. 3 7 1 7 \\ 0. 3 7 1 7 \\ 0. 6 0 1 5 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\lambda_ {3} \doteq 6. 8 5 4; \quad \phi_ {3} \doteq \left[ \begin{array}{c} - 0. 6 0 1 5 \\ 0. 3 7 1 7 \\ 0. 3 7 1 7 \\ - 0. 6 0 1 5 \end{array} \right]
|
|
$$
|
|
|
|
$$
|
|
\lambda_ {4} \doteq 1 3. 0 9; \quad \phi_ {4} \doteq \left[ \begin{array}{c} - 0. 3 7 1 7 \\ 0. 6 0 1 5 \\ - 0. 6 0 1 5 \\ 0. 3 7 1 7 \end{array} \right]
|
|
$$
|
|
|
|
It should be noted that the eigenvalues and eigenvectors did not appear in the usual order in the approximation for $\Lambda$ and $\Phi$ .
|
|
|
|
In the following example we demonstrate the quadratic convergence when the off-diagonal elements are already small (see J. H. Wilkinson [B]).
|
|
|
|
EXAMPLE 11.10: Consider the Jacobi solution of the eigenproblem $\mathbf{K}\phi = \lambda \phi$ , where
|
|
|
|
$$
|
|
\mathbf {K} = \left[ \begin{array}{l l l} k _ {1 1} & o (\epsilon) & o (\epsilon) \\ o (\epsilon) & k _ {2 2} & o (\epsilon) \\ o (\epsilon) & o (\epsilon) & k _ {3 3} \end{array} \right]
|
|
$$
|
|
|
|
The symbol $o(\epsilon)$ signifies "of order $\epsilon$ ," where $\epsilon \ll k_{ii}$ , $i = 1, 2, 3$ . Show that after one complete sweep, all off-diagonal elements are of order $\epsilon^2$ , meaning that convergence is quadratic.
|
|
|
|
Since the rotations to be applied are small, we make the assumption that $\sin \theta = \theta$ and $\cos \theta = 1$ . Hence, the relation in (11.80) gives
|
|
|
|
$$
|
|
\theta = \frac {k _ {i j} ^ {(k)}}{k _ {i i} ^ {(k)} - k _ {j j} ^ {(k)}}
|
|
$$
|
|
|
|
In one sweep we need to set to zero, in succession, all off-diagonal elements. Using $K_{1} = K$ , we obtain $K_{2}$ by zeroing element (1, 2) in $K_{1}$ ,
|
|
|
|
$$
|
|
\mathbf {K} _ {2} = \mathbf {P} _ {1} ^ {T} \mathbf {K} _ {1} \mathbf {P} _ {1}
|
|
$$
|
|
|
|
where $\mathbf{P}_{1}=\left[\begin{array}{ccc}1 & \frac{-o(\epsilon)}{k_{11}-k_{22}} & 0 \\ \frac{o(\epsilon)}{k_{11}-k_{22}} & 1 & 0 \\ 0 & 0 & 1\end{array}\right]$
|
|
|
|
<!-- source-page: 936 -->
|
|
|
|
Hence, $\mathbf{K}_2 = \begin{bmatrix} k_{11} + o(\epsilon^2) & 0 & o(\epsilon) \\ 0 & k_{22} + o(\epsilon^2) & o(\epsilon) \\ o(\epsilon) & o(\epsilon) & k_{33} \end{bmatrix}$
|
|
|
|
Similarly, we zero element (1, 3) in $\mathbf{K}_2$ to obtain $\mathbf{K}_3$ ,
|
|
|
|
$$
|
|
\mathbf {K} _ {3} = \left[ \begin{array}{c c c} k _ {1 1} + o \left(\epsilon^ {2}\right) & o \left(\epsilon^ {2}\right) & 0 \\ o \left(\epsilon^ {2}\right) & k _ {2 2} + o \left(\epsilon^ {2}\right) & o (\epsilon) \\ 0 & o (\epsilon) & k _ {3 3} + o \left(\epsilon^ {2}\right) \end{array} \right]
|
|
$$
|
|
|
|
Finally, we zero element (2, 3) in $\mathbf{K}_3$ and have
|
|
|
|
$$
|
|
\mathbf {K} _ {4} = \left[ \begin{array}{c c c} k _ {1 1} + o \left(\epsilon^ {2}\right) & o \left(\epsilon^ {2}\right) & o \left(\epsilon^ {2}\right) \\ o \left(\epsilon^ {2}\right) & k _ {2 2} + o \left(\epsilon^ {2}\right) & 0 \\ o \left(\epsilon^ {2}\right) & 0 & k _ {3 3} + o \left(\epsilon^ {2}\right) \end{array} \right]
|
|
$$
|
|
|
|
with all off-diagonal elements at least $o(\epsilon^{2})$ .
|
|
|
|
# 11.3.2 The Generalized Jacobi Method
|
|
|
|
In the previous section we discussed the solution of the standard eigenproblem $K\phi = \lambda\phi$ using the conventional Jacobi rotation matrices in order to reduce K to diagonal form. To solve the generalized problem $K\phi = \lambda M\phi$ , $M \neq I$ , using the standard Jacobi method, it would be necessary to first transform the problem into a standard form. However, this transformation can be dispensed with by using a generalized Jacobi solution method that operates directly on K and M (see S. Falk and P. Langemeyer [A] and K. J. Bathe [A]). The algorithm proceeds as summarized in (11.72) to (11.76) and is a natural extension of the standard Jacobi solution scheme; i.e., the generalized method reduces to the scheme presented for the problem $K\phi = \lambda\phi$ when M is an identity matrix.
|
|
|
|
Referring to the discussion in the previous section, in the generalized Jacobi iteration we use the following matrix $P_{k}$ :
|
|
|
|
$$
|
|
\mathbf {P} _ {k} = \left[ \begin{array}{c c c c c c c c c} 1 & & & & i \text {th} & & j \text {th column} \\ & \cdot & & & & & & & \\ & & \cdot & & & & & & \\ & & & \cdot & & & & & \\ & & & & 1 & & \alpha & & - i \text {th} \\ & & & & & \cdot & & & \\ & & & & & \gamma & & 1 & - j \text {th row} \\ & & & & & & & \cdot & \\ & & & & & & & \cdot & \\ & & & & & & & & 1 \end{array} \right] \tag {11.85}
|
|
$$
|
|
|
|
where the constants $\alpha$ and $\gamma$ are selected in such a way as to reduce to zero simultaneously elements $(i,j)$ in $\mathbf{K}_k$ and $\mathbf{M}_k$ . Therefore, the values of $\alpha$ and $\gamma$ are a function of the elements $k_{ij}^{(k)}, k_{ii}^{(k)}, k_{jj}^{(k)}, m_{ij}^{(k)}, m_{ii}^{(k)}$ , and $m_{jj}^{(k)}$ , where the superscript $(k)$ indicates that the $k$ th iteration is considered. Performing the multiplications $\mathbf{P}_k^T\mathbf{K}_k\mathbf{P}_k$ and $\mathbf{P}_k^T\mathbf{M}_k\mathbf{P}_k$ and using the condition
|
|
|
|
<!-- source-page: 937 -->
|
|
|
|
that $k_{ij}^{(k + 1)}$ and $m_{ij}^{(k + 1)}$ shall be zero, we obtain the following two equations for $\alpha$ and $\gamma$ :
|
|
|
|
$$
|
|
\alpha k _ {i i} ^ {(k)} + (1 + \alpha \gamma) k _ {i j} ^ {(k)} + \gamma k _ {j j} ^ {(k)} = 0 \tag {11.86}
|
|
$$
|
|
|
|
and $\alpha m_{ii}^{(k)} + (1 + \alpha \gamma)m_{ij}^{(k)} + \gamma m_{jj}^{(k)} = 0$ (11.87)
|
|
|
|
If $\frac{k_{ii}^{(k)}}{m_{ii}^{(k)}} = \frac{k_{jj}^{(k)}}{m_{jj}^{(k)}} = \frac{k_{ij}^{(k)}}{m_{ij}^{(k)}}$
|
|
|
|
(i.e., the submatrices considered are scalar multiples, which may be regarded to be a trivial case), we use $\alpha = 0$ and $\gamma = -k_{ij}^{(k)} / k_{jj}^{(k)}$ . In general, to solve for $\alpha$ and $\gamma$ from (11.86) and (11.87), we define
|
|
|
|
$$
|
|
\left. \begin{array}{l} \bar {k} _ {i i} ^ {(k)} = k _ {i i} ^ {(k)} m _ {i j} ^ {(k)} - m _ {i i} ^ {(k)} k _ {i j} ^ {(k)} \\ \bar {k} _ {j j} ^ {(k)} = k _ {j j} ^ {(k)} m _ {i j} ^ {(k)} - m _ {j j} ^ {(k)} k _ {i j} ^ {(k)} \\ \bar {k} ^ {(k)} = k _ {i i} ^ {(k)} m _ {j j} ^ {(k)} - k _ {j j} ^ {(k)} m _ {i i} ^ {(k)} \end{array} \right\} \tag {11.88}
|
|
$$
|
|
|
|
and $\gamma = -\frac{\bar{k}_{ii}^{(k)}}{x};\quad \alpha = \frac{\bar{k}_{jj}^{(k)}}{x} \tag{11.89}$
|
|
|
|
The value of x needed to obtain $\alpha$ and $\gamma$ is then to be determined using
|
|
|
|
$$
|
|
x = \frac {\bar {k} ^ {(k)}}{2} + \operatorname{sign} \left(\bar {k} ^ {(k)}\right) \sqrt {\left(\frac {\bar {k} ^ {(k)}}{2}\right) ^ {2} + \bar {k} _ {i i} ^ {(k)} \bar {k} _ {j j} ^ {(k)}} \tag {11.90}
|
|
$$
|
|
|
|
The relations for $\alpha$ and $\gamma$ are used and have primarily been developed for the case of M being a positive definite full or banded mass matrix. In that case (and, in fact, also under less restrictive conditions), we have
|
|
|
|
$$
|
|
\left(\frac {\bar {k} ^ {(k)}}{2}\right) ^ {2} + \bar {k} _ {i i} ^ {(k)} \bar {k} _ {j j} ^ {(k)} > 0
|
|
$$
|
|
|
|
and hence $x$ is always nonzero. In addition, $\det \mathbf{P}_k \neq 0$ , which indeed is the necessary condition for the algorithm to work.
|
|
|
|
The generalized Jacobi solution procedure has been used a great deal in the subspace iteration method (see Section 11.6) and when a consistent mass idealization is employed. However, other situations may arise as well. Assume that M is a diagonal mass matrix, $M \neq I$ and $m_{ii} > 0$ , in which case we employ in (11.88),
|
|
|
|
$$
|
|
\bar {k} _ {i i} ^ {(k)} = - m _ {i i} ^ {(k)} k _ {i j} ^ {(k)}; \quad \bar {k} _ {j j} ^ {(k)} = - m _ {j j} ^ {(k)} k _ {i j} ^ {(k)} \tag {11.91}
|
|
$$
|
|
|
|
and otherwise (11.85) to (11.90) are used as before. However, if $\mathbf{M} = \mathbf{I}$ , the relation in (11.87) yields $\alpha = -\gamma$ , and we recognize that $\mathbf{P}_k$ in (11.85) is a multiple of the rotation matrix defined in (11.79) (see Example 11.11). In addition, it should be mentioned that the solution procedure can be adapted to solve the problem $\mathbf{K}\phi = \lambda \mathbf{M}\phi$ when $\mathbf{M}$ is a diagonal matrix with some zero diagonal elements.
|
|
|
|
The complete solution process is analogous to the Jacobi iteration in the solution of the problem $K\phi = \lambda\phi$ , which was presented in the preceding section. The differences are that now a mass coupling factor $[(m_{ij}^{(k)})^{2}/m_{ii}^{(k)}m_{jj}^{(k)}]^{1/2}$ must also be calculated, unless M is diagonal, and the transformation is applied to $K_{k}$ and $M_{k}$ . Convergence is measured by comparing successive eigenvalue approximations and by testing if all off-diagonal elements are small enough; i.e., with l being the last iteration, convergence has been achieved if
|
|
|
|
$$
|
|
\frac {\left| \lambda_ {i} ^ {(l + 1)} - \lambda_ {i} ^ {(l)} \right|}{\lambda_ {i} ^ {(l + 1)}} \leq 1 0 ^ {- s}; \quad i = 1, \dots , n \tag {11.92}
|
|
$$
|
|
|
|
<!-- source-page: 938 -->
|
|
|
|
where
|
|
|
|
$$
|
|
\lambda_ {i} ^ {(l)} = \frac {k _ {i i} ^ {(l)}}{m _ {i i} ^ {(l)}}; \quad \lambda_ {i} ^ {(l + 1)} = \frac {k _ {i i} ^ {(l + 1)}}{m _ {i i} ^ {(l + 1)}} \tag {11.93}
|
|
$$
|
|
|
|
and
|
|
|
|
$$
|
|
\left[ \frac {(k _ {i j} ^ {(l + 1)}) ^ {2}}{k _ {i i} ^ {(l + 1)} k _ {j j} ^ {(l + 1)}} \right] ^ {1 / 2} \leq 1 0 ^ {- s}; \quad \left[ \frac {(m _ {i j} ^ {(l + 1)}) ^ {2}}{m _ {i i} ^ {(l + 1)} m _ {j j} ^ {(l + 1)}} \right] ^ {1 / 2} \leq 1 0 ^ {- s}; \quad \text { all } i, j; i < j \tag {11.94}
|
|
$$
|
|
|
|
where $10^{-s}$ is the convergence tolerance.
|
|
|
|
Table 11.1 summarizes the solution procedure for the case of M being full (or banded) and positive definite. The relations given in Table 11.1 are employed directly in the subroutine JACOBI, which is presented at the end of this section. Table 11.1 also gives an operation count of the solution process and the storage requirements. The total number of operations in one sweep as given in the table are an upper bound because it is assumed that both matrices are full and that all off-diagonal elements are zeroed; i.e., the threshold tolerance is never passed. Considering the number of sweeps required for solution, the same experience as with the solution of standard eigenproblems holds; i.e., with m = 2 and s = 12 in the iteration (see Section 11.3.1) about six sweeps are required for solution of the eigensystem to high accuracy.
|
|
|
|
TABLE 11.1 Summary of generalized Jacobi solution
|
|
<table><tr><td>Operation</td><td>Calculation</td><td>Number of operations</td><td>Required storage</td></tr><tr><td>Calculation of coupling factors</td><td> $\frac{(k_{ij}^{(k)})^2}{k_{ii}^{(k)} k_{jj}^{(k)}}; \frac{(m_{ij}^{(k)})^2}{m_{ii}^{(k)} m_{jj}^{(k)}}$ </td><td>6</td><td></td></tr><tr><td rowspan="3">Transformation to zero elements (i, j)</td><td> $\overline{k}_{ii}^{(k)} = k_{ii}^{(k)} m_{ij}^{(k)} - m_{ii}^{(k)} k_{ij}^{(k)}$ </td><td></td><td></td></tr><tr><td> $\overline{k}_{jj}^{(k)} = k_{jj}^{(k)} m_{ij}^{(k)} - m_{jj}^{(k)} k_{ij}^{(k)}$ </td><td></td><td></td></tr><tr><td> $\overline{k}^{(k)} = k_{ii}^{(k)} m_{jj}^{(k)} - k_{jj}^{(k)} m_{ii}^{(k)}$ </td><td></td><td></td></tr><tr><td rowspan="3"></td><td> $x = \frac{\overline{k}^{(k)}}{2} + (\text{sign } \overline{k}^{(k)}) \sqrt{\left( \frac{\overline{k}^{(k)}}{2} \right)^2 + \overline{k}_{ii}^{(k)} \overline{k}_{jj}^{(k)}}$ </td><td>4n + 12</td><td>Using symmetry of matrices $n(n + 2)$ </td></tr><tr><td> $\gamma = -\frac{\overline{k}_{ii}^{(k)}}{x}, \alpha = \frac{\overline{k}_{jj}^{(k)}}{x}$ </td><td></td><td></td></tr><tr><td> $\mathbf{K}_{k+1} = \mathbf{P}_k^T \mathbf{K}_k \mathbf{P}_k, \quad \mathbf{M}_{k+1} = \mathbf{P}_k^T \mathbf{M}_k \mathbf{P}_k$ </td><td></td><td></td></tr><tr><td>Calculation of eigenvectors</td><td> $(\mathbf{P}_1 \ldots \mathbf{P}_{k-1}) \mathbf{P}_k$ </td><td>2n</td><td> $n^2$ </td></tr><tr><td>Total for one sweep</td><td></td><td> $3n^3 + 6n^2$ </td><td> $2n^2 + 2n$ </td></tr></table>
|
|
|
|
The following examples demonstrate some of the characteristics of the generalized Jacobi solution algorithm.
|
|
|
|
EXAMPLE 11.11: Prove that the generalized Jacobi method reduces to the standard technique when $\mathbf{M} = \mathbf{I}$ .
|
|
|
|
For the proof we need only consider the calculation of the transformation matrices that would be used to zero typical off-diagonal elements. We want to show that the transformation matrices obtained in the standard and generalized Jacobi methods are multiples of each other; namely, in that case we could, by proper scaling, obtain the standard method from the
|
|
|
|
<!-- source-page: 939 -->
|
|
|
|
generalized scheme. Since each step of iteration consists of applying a rotation in the $(i,j)$ th plane, we can without loss of generality consider the solution of the problem
|
|
|
|
$$
|
|
\left[ \begin{array}{l l} k _ {1 1} & k _ {1 2} \\ k _ {1 2} & k _ {2 2} \end{array} \right] \boldsymbol {\phi} = \lambda \boldsymbol {\phi}
|
|
$$
|
|
|
|
Using (11.88) to (11.90), we thus obtain
|
|
|
|
$$
|
|
\alpha = - \gamma ; \quad \mathbf {P} _ {1} = \left[ \begin{array}{c c} 1 & - \gamma \\ \gamma & 1 \end{array} \right] \tag {a}
|
|
$$
|
|
|
|
and $\gamma = \frac{-k_{11} + k_{22} \pm \sqrt{(k_{11} - k_{22})^2 + 4k_{12}^2}}{2k_{12}}$
|
|
|
|
On the other hand, in the standard Jacobi solution we use
|
|
|
|
$$
|
|
\mathbf {P} _ {1} = \left[ \begin{array}{c c} \cos \theta & - \sin \theta \\ \sin \theta & \cos \theta \end{array} \right]
|
|
$$
|
|
|
|
which may be written as
|
|
|
|
$$
|
|
\mathbf {P} _ {1} = \cos \theta \left[ \begin{array}{c c} 1 & - \tan \theta \\ \tan \theta & 1 \end{array} \right] \tag {b}
|
|
$$
|
|
|
|
Thus, $\mathbf{P}_1$ in (b) would be a multiple of $\mathbf{P}_1$ in (a) if $\tan \theta = \gamma$ . In the standard Jacobi method we obtain $\tan 2\theta$ using (11.80). In this case, we have
|
|
|
|
$$
|
|
\tan 2 \theta = \frac {2 k _ {1 2}}{k _ {1 1} - k _ {2 2}} \tag {c}
|
|
$$
|
|
|
|
We also have, by simple trigonometry,
|
|
|
|
$$
|
|
\tan 2 \theta = \frac {2 \tan \theta}{1 - \tan^ {2} \theta} \tag {d}
|
|
$$
|
|
|
|
Using (c) and (d), we can solve for $\tan \theta$ to be used in (b) and obtain
|
|
|
|
$$
|
|
\tan \theta = \frac {- k _ {1 1} + k _ {2 2} \pm \sqrt {\left(k _ {1 1} - k _ {2 2}\right) ^ {2} + 4 k _ {1 2} ^ {2}}}{2 k _ {1 2}}
|
|
$$
|
|
|
|
Hence, $\gamma = \tan \theta$ and the generalized Jacobi iteration is equivalent to the standard method when $\mathbf{M} = \mathbf{I}$ .
|
|
|
|
EXAMPLE 11.12: Use the generalized Jacobi method to calculate the eigensystem of the problem $\mathbf{K}\phi = \lambda \mathbf{M}\phi$ .
|
|
|
|
(1) In the first case let
|
|
|
|
$$
|
|
\mathbf {K} = \left[ \begin{array}{c c} 1 & - 1 \\ - 1 & 1 \end{array} \right]; \quad \mathbf {M} = \left[ \begin{array}{c c} 2 & 1 \\ 1 & 2 \end{array} \right]
|
|
$$
|
|
|
|
We note that $\mathbf{K}$ is singular, and hence we expect a zero eigenvalue.
|
|
|
|
(2) Then let
|
|
|
|
$$
|
|
\mathbf {K} = \left[ \begin{array}{l l} 2 & 1 \\ 1 & 2 \end{array} \right]; \quad \mathbf {M} = \left[ \begin{array}{l l} 2 & 0 \\ 0 & 0 \end{array} \right]
|
|
$$
|
|
|
|
in which case we have an infinite eigenvalue.
|
|
|
|
<!-- source-page: 940 -->
|
|
|
|
For solution, we use the relations in (11.85) to (11.90). Considering the problem in case (1), we obtain
|
|
|
|
$$
|
|
\begin{array}{l} \bar {k} _ {1 1} ^ {(1)} = 3; \quad \bar {k} _ {2 2} ^ {(1)} = 3; \quad \bar {k} ^ {(1)} = 0 \\ x = 3; \quad \gamma = - 1; \quad \alpha = 1 \\ \mathbf {P} _ {1} = \left[ \begin{array}{c c} 1 & 1 \\ - 1 & 1 \end{array} \right] \\ \end{array}
|
|
$$
|
|
|
|
Hence, $\mathbf{P}_1^T\mathbf{K}\mathbf{P}_1 = \begin{bmatrix} 4 & 0 \\ 0 & 0 \end{bmatrix}; \quad \mathbf{P}_1^T\mathbf{M}\mathbf{P}_1 = \begin{bmatrix} 2 & 0 \\ 0 & 6 \end{bmatrix}$
|
|
|
|
To obtain $\Lambda$ and $\Phi$ we use (11.75) and (11.76) and arrange the columns in the matrices in the appropriate order. Hence,
|
|
|
|
$$
|
|
\boldsymbol {\Lambda} = \left[ \begin{array}{l l} 0 & \\ & 2 \end{array} \right]; \quad \boldsymbol {\Phi} = \left[ \begin{array}{c c} \frac {1}{\sqrt {6}} & \frac {1}{\sqrt {2}} \\ \frac {1}{\sqrt {6}} & - \frac {1}{\sqrt {2}} \end{array} \right]
|
|
$$
|
|
|
|
Now consider the problem in case (2). Here we have
|
|
|
|
$$
|
|
\begin{array}{l} \bar {k} _ {1 1} ^ {(1)} = - 2; \quad \bar {k} _ {2 2} ^ {(1)} = 0; \quad \bar {k} ^ {(1)} = - 4 \\ x = - 4; \quad \alpha = 0; \quad \gamma = - \frac {1}{2} \\ \mathbf {P} _ {1} = \left[ \begin{array}{c c} 1 & 0 \\ - \frac {1}{2} & 1 \end{array} \right] \\ \end{array}
|
|
$$
|
|
|
|
Hence, $\mathbf{P}_1^T\mathbf{K}\mathbf{P}_1 = \begin{bmatrix} \frac{3}{2} & 0\\ 0 & 2 \end{bmatrix};\qquad \mathbf{P}_1^T\mathbf{M}\mathbf{P}_1 = \begin{bmatrix} 2 & 0\\ 0 & 0 \end{bmatrix}$
|
|
|
|
and $\Lambda = \begin{bmatrix} \frac{3}{4} & \\ & \infty \end{bmatrix};\quad \phi_{1} = \begin{bmatrix} \frac{1}{\sqrt{2}} \\ -\frac{1}{2\sqrt{2}} \end{bmatrix}$
|
|
|
|
The above discussion of the generalized Jacobi solution method has already indicated in some way the advantages of the solution technique. First, the transformation of the generalized eigenproblem to the standard form is avoided. This is particularly advantageous (1) when the matrices are ill-conditioned, and (2) when the off-diagonal elements in K and M are already small or, equivalently when there are only a few nonzero off-diagonal elements. In the first case the direct solution of $K\phi = \lambda M\phi$ avoids the solution of a standard eigenproblem of a matrix with very large and very small elements (see Section 10.2.5). In the second case the eigenproblem is already nearly solved, because the zeroing of small or only a few off-diagonal elements in K and M will not result in a large change in the diagonal elements of the matrices, the ratios of which are the eigenvalues. In addition, fast convergence can be expected when the off-diagonal elements are small (see Section 11.3.1). We will see that this case arises in the subspace iteration method described in Section 11.6, which is one reason why the generalized Jacobi method is used effectively in that technique.
|