CFD Online URL
[Sponsors]
Home > Wiki > Gauss-Seidel method

Gauss-Seidel method

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
(Added an example calculation (finally))
(Example Calculation: Added more data)
Line 38: Line 38:
In some cases, we need not even explicitly represent the matrix we are solving.  Consider the simple heat equation problem  
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>
+
:<math>\nabla^2 T(x) = 0,\ x\in [0,1]</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
+
subject to the boundary conditions <math>T(0)=0</math> and <math>T(1)=1</math>.  The exact solution to this problem is <math>T(x)=x</math>.  The standard second-order finite difference discretization is
:<math> T_{i-1}-2T_i+T_{i+1} = 0,</math>
:<math> T_{i-1}-2T_i+T_{i+1} = 0,</math>
Line 84: Line 84:
:<math> T_i = \frac{1}{2}(T_{i-1}+T_{i+1}).</math>
:<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.
+
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 <math>L_2</math> norm error.
   
   
{| align=center border=1
{| align=center border=1
-
! Iteration !! T(1) !! T(2) !! T(3) !! T(4) !! T(5) !! L2 error  
+
|+ Gauss-Seidel Solution (i=2,3,4)
 +
! Iteration !! <math>T_1</math> !! <math>T_2</math> !! <math>T_3</math> !! <math>T_4</math> !! <math>T_5</math> !! <math>L_2</math> error  
|-
|-
|    0  
|    0  
Line 177: Line 178:
|  1.4648E-03
|  1.4648E-03
|}
|}
 +
 +
In this situation, the direction that we "sweep" is important - if we step though the solution vector in the opposite direction, the solution moves away from the chosen initial condition (zero everywhere in the interior) more quickly.  This gives us (slightly) faster convergence, as shown in the table below.
 +
 +
{| align=center border=1
 +
|+ Gauss-Seidel Solution (i=4,3,2)
 +
! Iteration !! <math>T_1</math> !! <math>T_2</math> !! <math>T_3</math> !! <math>T_4</math> !! <math>T_5</math> !! <math>L_2</math> error
 +
|-
 +
|    0
 +
|  0.0000E+00
 +
|  0.0000E+00
 +
|  0.0000E+00
 +
|  0.0000E+00
 +
|  1.0000E+00
 +
|  1.0000E+00
 +
|-
 +
|    1
 +
|  0.0000E+00
 +
|  1.2500E-01
 +
|  2.5000E-01
 +
|  5.0000E-01
 +
|  1.0000E+00
 +
|  3.7500E-01
 +
|-
 +
|    2
 +
|  0.0000E+00
 +
|  1.8750E-01
 +
|  3.7500E-01
 +
|  6.2500E-01
 +
|  1.0000E+00
 +
|  1.8750E-01
 +
|-
 +
|    3
 +
|  0.0000E+00
 +
|  2.1875E-01
 +
|  4.3750E-01
 +
|  6.8750E-01
 +
|  1.0000E+00
 +
|  9.3750E-02
 +
|-
 +
|    4
 +
|  0.0000E+00
 +
|  2.3438E-01
 +
|  4.6875E-01
 +
|  7.1875E-01
 +
|  1.0000E+00
 +
|  4.6875E-02
 +
|-
 +
|    5
 +
|  0.0000E+00
 +
|  2.4219E-01
 +
|  4.8438E-01
 +
|  7.3438E-01
 +
|  1.0000E+00
 +
|  2.3438E-02
 +
|-
 +
|    6
 +
|  0.0000E+00
 +
|  2.4609E-01
 +
|  4.9219E-01
 +
|  7.4219E-01
 +
|  1.0000E+00
 +
|  1.1719E-02
 +
|-
 +
|    7
 +
|  0.0000E+00
 +
|  2.4805E-01
 +
|  4.9609E-01
 +
|  7.4609E-01
 +
|  1.0000E+00
 +
|  5.8594E-03
 +
|-
 +
|    8
 +
|  0.0000E+00
 +
|  2.4902E-01
 +
|  4.9805E-01
 +
|  7.4805E-01
 +
|  1.0000E+00
 +
|  2.9297E-03
 +
|-
 +
|    9
 +
|  0.0000E+00
 +
|  2.4951E-01
 +
|  4.9902E-01
 +
|  7.4902E-01
 +
|  1.0000E+00
 +
|  1.4648E-03
 +
|-
 +
|  10
 +
|  0.0000E+00
 +
|  2.4976E-01
 +
|  4.9951E-01
 +
|  7.4951E-01
 +
|  1.0000E+00
 +
|  7.3242E-04
 +
|}
 +
 +
For this toy example, there is not large penalty for choosing the wrong sweep direction.  For some of the more complicated variants of Gauss-Seidel, there is a substantial penalty - the sweep direction determines (in a vague sense) the direction in which information travels.

Revision as of 02:45, 8 April 2006

Introduction

We seek the solution to set of linear equations:

 A  \phi = b

In matrix terms, the the Gauss-Seidel iteration can be expressed as

 
\phi^{(k+1)}  = \left( {D - L} \right)^{ - 1} \left( {U\phi^{(k)}  + b} \right),

where D, L, and U represent the diagonal, lower triangular, and upper triangular parts of the coefficient matrix A and k 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:

 
\phi^{(k+1)}_i  = \frac{1}{a_{ii}} \left(b_i - \sum_{j<i}a_{ij}\phi^{(k+1)}_j-\sum_{j>i}a_{ij}\phi^{(k)}_j\right),\, i=1,2,\ldots,n.

Note that the computation of \phi^{(k+1)}_i uses only those elements of \phi^{(k+1)} that have already been computed and only those elements of \phi^{(k)} that have yet to be advanced to iteration k+1. This means that no additional storage is required, and the computation can be done in place (\phi^{(k+1)} replaces \phi^{(k)}). 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 \phi^{0}
for k := 1 step 1 until convergence do
for i := 1 step until n do
 \sigma = 0
for j := 1 step until i-1 do
 \sigma  = \sigma  + a_{ij} \phi_j^{(k)}
end (j-loop)
for j := i+1 step until n do
 \sigma  = \sigma  + a_{ij} \phi_j^{(k-1)}
end (j-loop)
  \phi_i^{(k)}  = {{\left( {b_i  - \sigma } \right)} \over {a_{ii} }}
end (i-loop)
check if convergence is reached
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

\nabla^2 T(x) = 0,\ x\in [0,1]

subject to the boundary conditions T(0)=0 and T(1)=1. The exact solution to this problem is T(x)=x. The standard second-order finite difference discretization is

 T_{i-1}-2T_i+T_{i+1} = 0,

where T_i is the (discrete) solution available at uniformly spaced nodes. In matrix terms, this can be written as

 
\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].

However, for any given T_i for 1 < i < n, we can write

 T_i = \frac{1}{2}(T_{i-1}+T_{i+1}).

Then, stepping through the solution vector from i=2 to i=n-1, we can compute the next iterate from the two surrounding values. Note that (in this scheme), T_{i+1} is from the previous iteration, while T_{i-1} 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 L_2 norm error.

Gauss-Seidel Solution (i=2,3,4)
Iteration T_1 T_2 T_3 T_4 T_5 L_2 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

In this situation, the direction that we "sweep" is important - if we step though the solution vector in the opposite direction, the solution moves away from the chosen initial condition (zero everywhere in the interior) more quickly. This gives us (slightly) faster convergence, as shown in the table below.

Gauss-Seidel Solution (i=4,3,2)
Iteration T_1 T_2 T_3 T_4 T_5 L_2 error
0 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 1.0000E+00 1.0000E+00
1 0.0000E+00 1.2500E-01 2.5000E-01 5.0000E-01 1.0000E+00 3.7500E-01
2 0.0000E+00 1.8750E-01 3.7500E-01 6.2500E-01 1.0000E+00 1.8750E-01
3 0.0000E+00 2.1875E-01 4.3750E-01 6.8750E-01 1.0000E+00 9.3750E-02
4 0.0000E+00 2.3438E-01 4.6875E-01 7.1875E-01 1.0000E+00 4.6875E-02
5 0.0000E+00 2.4219E-01 4.8438E-01 7.3438E-01 1.0000E+00 2.3438E-02
6 0.0000E+00 2.4609E-01 4.9219E-01 7.4219E-01 1.0000E+00 1.1719E-02
7 0.0000E+00 2.4805E-01 4.9609E-01 7.4609E-01 1.0000E+00 5.8594E-03
8 0.0000E+00 2.4902E-01 4.9805E-01 7.4805E-01 1.0000E+00 2.9297E-03
9 0.0000E+00 2.4951E-01 4.9902E-01 7.4902E-01 1.0000E+00 1.4648E-03
10 0.0000E+00 2.4976E-01 4.9951E-01 7.4951E-01 1.0000E+00 7.3242E-04

For this toy example, there is not large penalty for choosing the wrong sweep direction. For some of the more complicated variants of Gauss-Seidel, there is a substantial penalty - the sweep direction determines (in a vague sense) the direction in which information travels.

My wiki