Gauss-Seidel method
From CFD-Wiki
m (Removed Wikipedia link - that article was based on this one.) |
(Added an example calculation (finally)) |
||
Line 34: | Line 34: | ||
:: check if convergence is reached | :: check if convergence is reached | ||
: end (k-loop) | : end (k-loop) | ||
+ | |||
+ | == Example Calculation == | ||
+ | In some cases, we need not even explicitly represent the matrix we are solving. Consider the simple heat equation problem | ||
+ | |||
+ | :<math>\nabla^2 T = 0</math> | ||
+ | |||
+ | on the unit interval subject to the boundary conditions <math>T(0)=0</math> and <math>T(1)=1</math>. The standard second-order finite difference discretization is | ||
+ | |||
+ | :<math> T_{i-1}-2T_i+T_{i+1} = 0,</math> | ||
+ | |||
+ | where <math>T_i</math> is the (discrete) solution available at uniformly spaced nodes. In matrix terms, this can be written as | ||
+ | |||
+ | :<math> | ||
+ | \left[ | ||
+ | \begin{matrix} | ||
+ | {1 } & {0 } & {0 } & \cdot & \cdot & { 0 } \\ | ||
+ | {1 } & {-2 } & {1 } & { } & { } & \cdot \\ | ||
+ | { 0 } & {1 } & {-2 } & { 1 } & { } & \cdot \\ | ||
+ | { } & { } & \cdot & \cdot & \cdot & \cdot \\ | ||
+ | { } & { } & { } & {1 } & {-2 } & {1 }\\ | ||
+ | { 0 } & \cdot & \cdot & { 0 } & { 0 } & {1 }\\ | ||
+ | \end{matrix} | ||
+ | \right] | ||
+ | \left[ | ||
+ | \begin{matrix} | ||
+ | {T_1 } \\ | ||
+ | {T_2 } \\ | ||
+ | {T_3 } \\ | ||
+ | \cdot \\ | ||
+ | {T_{n-1} } \\ | ||
+ | {T_n } \\ | ||
+ | \end{matrix} | ||
+ | \right] | ||
+ | = | ||
+ | \left[ | ||
+ | \begin{matrix} | ||
+ | {0} \\ | ||
+ | {0} \\ | ||
+ | {0} \\ | ||
+ | \cdot \\ | ||
+ | {0} \\ | ||
+ | {1} \\ | ||
+ | \end{matrix} | ||
+ | \right]. | ||
+ | </math> | ||
+ | |||
+ | However, for any given <math>T_i</math> for <math>1 < i < n</math>, we can write | ||
+ | |||
+ | :<math> T_i = \frac{1}{2}(T_{i-1}+T_{i+1}).</math> | ||
+ | |||
+ | Then, stepping through the solution vector from <math>i=2</math> to <math>i=n-1</math>, we can compute the next iterate from the two surrounding values. Note that (in this scheme), <math>T_{i+1}</math> is from the previous iteration, while <math>T_{i-1}</math> is from the current iteration. The following table gives the results of 10 iterations with 5 nodes (3 interior and 2 boundary) as well as L2 norm error. | ||
+ | |||
+ | {| align=center border=1 | ||
+ | ! Iteration !! T(1) !! T(2) !! T(3) !! T(4) !! T(5) !! L2 error | ||
+ | |- | ||
+ | | 0 | ||
+ | | 0.0000E+00 | ||
+ | | 0.0000E+00 | ||
+ | | 0.0000E+00 | ||
+ | | 0.0000E+00 | ||
+ | | 1.0000E+00 | ||
+ | | 1.0000E+00 | ||
+ | |- | ||
+ | | 1 | ||
+ | | 0.0000E+00 | ||
+ | | 0.0000E+00 | ||
+ | | 0.0000E+00 | ||
+ | | 5.0000E-01 | ||
+ | | 1.0000E+00 | ||
+ | | 6.1237E-01 | ||
+ | |- | ||
+ | | 2 | ||
+ | | 0.0000E+00 | ||
+ | | 0.0000E+00 | ||
+ | | 2.5000E-01 | ||
+ | | 6.2500E-01 | ||
+ | | 1.0000E+00 | ||
+ | | 3.7500E-01 | ||
+ | |- | ||
+ | | 3 | ||
+ | | 0.0000E+00 | ||
+ | | 1.2500E-01 | ||
+ | | 3.7500E-01 | ||
+ | | 6.8750E-01 | ||
+ | | 1.0000E+00 | ||
+ | | 1.8750E-01 | ||
+ | |- | ||
+ | | 4 | ||
+ | | 0.0000E+00 | ||
+ | | 1.8750E-01 | ||
+ | | 4.3750E-01 | ||
+ | | 7.1875E-01 | ||
+ | | 1.0000E+00 | ||
+ | | 9.3750E-02 | ||
+ | |- | ||
+ | | 5 | ||
+ | | 0.0000E+00 | ||
+ | | 2.1875E-01 | ||
+ | | 4.6875E-01 | ||
+ | | 7.3438E-01 | ||
+ | | 1.0000E+00 | ||
+ | | 4.6875E-02 | ||
+ | |- | ||
+ | | 6 | ||
+ | | 0.0000E+00 | ||
+ | | 2.3438E-01 | ||
+ | | 4.8438E-01 | ||
+ | | 7.4219E-01 | ||
+ | | 1.0000E+00 | ||
+ | | 2.3438E-02 | ||
+ | |- | ||
+ | | 7 | ||
+ | | 0.0000E+00 | ||
+ | | 2.4219E-01 | ||
+ | | 4.9219E-01 | ||
+ | | 7.4609E-01 | ||
+ | | 1.0000E+00 | ||
+ | | 1.1719E-02 | ||
+ | |- | ||
+ | | 8 | ||
+ | | 0.0000E+00 | ||
+ | | 2.4609E-01 | ||
+ | | 4.9609E-01 | ||
+ | | 7.4805E-01 | ||
+ | | 1.0000E+00 | ||
+ | | 5.8594E-03 | ||
+ | |- | ||
+ | | 9 | ||
+ | | 0.0000E+00 | ||
+ | | 2.4805E-01 | ||
+ | | 4.9805E-01 | ||
+ | | 7.4902E-01 | ||
+ | | 1.0000E+00 | ||
+ | | 2.9297E-03 | ||
+ | |- | ||
+ | | 10 | ||
+ | | 0.0000E+00 | ||
+ | | 2.4902E-01 | ||
+ | | 4.9902E-01 | ||
+ | | 7.4951E-01 | ||
+ | | 1.0000E+00 | ||
+ | | 1.4648E-03 | ||
+ | |} |
Revision as of 23:36, 7 April 2006
Introduction
We seek the solution to set of linear equations:
In matrix terms, the the Gauss-Seidel iteration can be expressed as
where , , and represent the diagonal, lower triangular, and upper triangular parts of the coefficient matrix and is the iteration count. This matrix expression is mainly of academic interest, and is not used to program the method. Rather, an element-based approach is used:
Note that the computation of uses only those elements of that have already been computed and only those elements of that have yet to be advanced to iteration . This means that no additional storage is required, and the computation can be done in place ( replaces ). While this might seem like a rather minor concern, for large systems it is unlikely that every iteration can be stored. Thus, unlike the Jacobi method, we do not have to do any vector copying should we wish to use only one storage vector. The iteration is generally continued until the changes made by an iteration are below some tolerance.
Algorithm
- Chose an initial guess
- for k := 1 step 1 until convergence do
- for i := 1 step until n do
-
- for j := 1 step until i-1 do
- end (j-loop)
- for j := i+1 step until n do
- end (j-loop)
-
- end (i-loop)
- check if convergence is reached
- for i := 1 step until n do
- end (k-loop)
Example Calculation
In some cases, we need not even explicitly represent the matrix we are solving. Consider the simple heat equation problem
on the unit interval subject to the boundary conditions and . The standard second-order finite difference discretization is
where is the (discrete) solution available at uniformly spaced nodes. In matrix terms, this can be written as
However, for any given for , we can write
Then, stepping through the solution vector from to , we can compute the next iterate from the two surrounding values. Note that (in this scheme), is from the previous iteration, while is from the current iteration. The following table gives the results of 10 iterations with 5 nodes (3 interior and 2 boundary) as well as L2 norm error.
Iteration | T(1) | T(2) | T(3) | T(4) | T(5) | L2 error |
---|---|---|---|---|---|---|
0 | 0.0000E+00 | 0.0000E+00 | 0.0000E+00 | 0.0000E+00 | 1.0000E+00 | 1.0000E+00 |
1 | 0.0000E+00 | 0.0000E+00 | 0.0000E+00 | 5.0000E-01 | 1.0000E+00 | 6.1237E-01 |
2 | 0.0000E+00 | 0.0000E+00 | 2.5000E-01 | 6.2500E-01 | 1.0000E+00 | 3.7500E-01 |
3 | 0.0000E+00 | 1.2500E-01 | 3.7500E-01 | 6.8750E-01 | 1.0000E+00 | 1.8750E-01 |
4 | 0.0000E+00 | 1.8750E-01 | 4.3750E-01 | 7.1875E-01 | 1.0000E+00 | 9.3750E-02 |
5 | 0.0000E+00 | 2.1875E-01 | 4.6875E-01 | 7.3438E-01 | 1.0000E+00 | 4.6875E-02 |
6 | 0.0000E+00 | 2.3438E-01 | 4.8438E-01 | 7.4219E-01 | 1.0000E+00 | 2.3438E-02 |
7 | 0.0000E+00 | 2.4219E-01 | 4.9219E-01 | 7.4609E-01 | 1.0000E+00 | 1.1719E-02 |
8 | 0.0000E+00 | 2.4609E-01 | 4.9609E-01 | 7.4805E-01 | 1.0000E+00 | 5.8594E-03 |
9 | 0.0000E+00 | 2.4805E-01 | 4.9805E-01 | 7.4902E-01 | 1.0000E+00 | 2.9297E-03 |
10 | 0.0000E+00 | 2.4902E-01 | 4.9902E-01 | 7.4951E-01 | 1.0000E+00 | 1.4648E-03 |