I know that for tridiagonal matrices the two iterative methods for linear system solving, the Gauss-Seidel method and the Jacobi one, either both converge or neither converges, and the Gauss-Seidel method converges twice as fast as the Jacobi one. Are there theorems relating the convergence speeds of these two methods when the system's matrix is not tridiagonal?
[Math] Jacobi vs. Gauss-Seidel: convergence
linear algebranumerical linear algebra
Related Solutions
I have found solution for my question and implemented it in C++:
/*
* www.twitter.com/torayeff
*
* Gauss-Seidel method:
* http://math.semestr.ru/optim/methodzeidel.php
* http://en.wikipedia.org/wiki/Gauss–Seidel_method
* http://ru.wikipedia.org/wiki/Метод_Гаусса_—_Зейделя
*
* To solve A*x = b transform to
* A(t)*A*x = A(t)*b <---> C*x = d
* (A(t) transpose of matrix A)
* (in this state, Gauss-Seidel method always congerges)
*/
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
#define MAXN 100
int N; //matrix dimension
double eps = 0.001;
double err = eps;
double A[MAXN][MAXN];
double b[MAXN][MAXN];
double At[MAXN][MAXN];
double C[MAXN][MAXN];
double d[MAXN][MAXN];
//convergence check
bool converge(double *xk, double *xkp)
{
double norm = 0;
for(int i=0; i<N; i++)
{
norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]);
}
if(sqrt(norm)>=eps)
return false;
err = eps;
return true;
}
int main()
{
printf("Input dimension N: ");
scanf("%d", &N);
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
double tmp;
printf("A[%d][%d] = ", i+1, j+1);
scanf("%lf", &tmp);
A[i][j] = tmp;
}
}
for(int i=0; i<N; i++)
{
double tmp;
printf("b[%d] = ", i+1);
scanf("%lf", &tmp);
b[i][0] = tmp;
}
//transpose A
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
At[i][j] = A[j][i];
}
}
//multiply At*A = C
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
C[i][j] = 0;
for(int k=0; k<N; k++)
{
C[i][j] = C[i][j] + At[i][k]*A[k][j];
}
}
}
//multiply At*b = d
for(int i=0; i<N; i++)
{
for(int j=0; j<1; j++)
{
d[i][j] = 0;
for(int k=0; k<N; k++)
{
d[i][j] = d[i][j] + At[i][k]*b[k][j];
}
}
}
//Solve Gauss-Seidel method for C*x = d
double x[N]; //current values
double p[N]; //previous values
for(int i=0; i<N; i++)
{
x[i] = 0;
p[i] = 0;
}
do
{
for(int i = 0; i<N; i++)
p[i] = x[i];
for(int i=0; i<N; i++)
{
double v = 0.0;
for(int j=0; j<i; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*x[j];
}
for(int j=i; j<N; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*p[j];
}
v += 1.0*d[i][0]/C[i][i];
x[i] = v;
}
}while (!converge(x, p));
//print solution
printf("Err. val.: %lf\n", err);
for(int i=0; i<N; i++)
{
printf("%lf ", x[i]);
}
return 0;
}
- $$ \begin{pmatrix} 1 & 0 & 1 \\ -1 & 1 & 0 \\ 1 & 2 & -3 \\ \end{pmatrix} $$
is convergent by J but divergent by GS. and $$ \begin{pmatrix} 1 & 0.5 & 0.5 \\ 0.5 & 1 & 0.5 \\ 0.5 & 0.5 & 1\\ \end{pmatrix} $$ is convergent by GS but divergent by J.
- In general,the convergent speed by GS is faster than J if them are also convergent.
Best Answer
Yes. From pg. 233 of A Friendly Introduction to Numerical Analysis by Brian Bradie, we have the following:
Theorem. Suppose $A$ is a an $n \times n$ matrix. If $ a_{ii} > 0$ for each $i$ and $a_{ij} \leq 0$ whenever $ i \neq j$, then one and only one of the following statements holds:
$0 \leq \rho (T_{gs}) < \rho (T_{jac}) < 1$
$1 < \rho (T_{jac}) < \rho (T_{gs})$
$\rho (T_{jac}) = \rho (T_{gs}) = 0$
$\rho (T_{jac}) =\rho (T_{gs}) = 1$
where $T$ is the iteration matrix that arises for each method, and $\rho (T)$ denotes the spectral radius of $T$.
Thus, for this larger class of matrices, the methods converge and diverge together. When they converge, Gauss-Seidel converges faster; when they diverge, Gauss-Seidel again does so faster.
If you want the proof of this, Bradie cites the following sources:
A. Ralston and P. Rabinowitz, A First Course in Numerical Analysis, 2nd edition, McGraw-Hill, New York, 1978.
D. M. Young, Iterative Solution of Large Linear Systems, Academic Press, New York, 1971.