SOR和SSOR迭代

计算数学与科学工程计算研究所 陆嵩

基本思想

对于一般的大型稀疏方程组的求解，我们一般使用迭代方法求解，而古典的迭代方法中，松弛方法最优。对于Jacobi迭代，Gauss-Seidel迭代和（对称）超松弛方法迭代，我们可以简单地用矩阵分裂的思想来考虑这个问题，也可以从方程组的角度来理解这个问题。Jacobi迭代是将第i个方程中的第i个分量单独挪到一边，每个方程单独迭代的过程。Gauss-Seidel迭代是在Jacobi迭代过程中，对于迭代变量，求解方程中若某分量产生了新值用新值而不用旧值的过程，SOR迭代是在新值和旧值之间取个加权平均的过程。SSOR呢，是对于方程组，先用SOR方法从前往后（下三角用新值）算一遍，再用SOR方法从后往前（上三角用新值）算一遍的当成一次迭代的过程。从方程推导出的矩阵和分裂方法得到的结果是一样的。SSOR比起SOR更有可能收敛，并且对于权值来说不是太敏感。

一般来说SOR（包含了Gauss-Seidel）方法是优于Jacobi方法的，但是由于Jacabi良好的并行性质，使其仍存于世上。

对于对称正定矩阵，Jacobi迭代在2D−A正定时收敛，Gauss-Seidel方法收敛，SOR和SSOR收敛（0<ω<2）。
对于严格对角占优和不可分弱严格对角占优矩阵（对角线不为零的是H矩阵），Jacobi迭代和Gauss-Seidel迭代收敛，SOR和SSOR在0<ω<2/[1+ρ(|J|)]时收敛。

让人比较容易记住的是，矩阵分裂方法中，Jacobi迭代矩阵分裂出的M为对角矩阵，Gauss-Seidel迭代矩阵分裂出的M为下三角，SOR方法的迭代矩阵分裂出的M为下三角在对角元上再乘上1/ω，与其对称的一个分裂出的M是上三角且对角元素乘以1/ω。

需要指出的是，虽然SSOR比起SOR更稳定，收敛的可能性更大，但是有时候，其实在SOR收敛的情况下，最优SOR收敛可能比最优SSOR收敛得要快。所以，在实际使用中，我们还应该根据实际情况来选择迭代方式。对于非对称正定，对对角占优的矩阵，这些方法一般不收敛。

程序与结果

• C代码

  #include<stdlib.h>#include<stdio.h>#include<math.h>#include<string.h>#define N 5#define w 1#define epsilon 0.0000001void main(){    double normD(double x0[N], double x1[N]);    void SORiteration1(double A[N][N], double b[N], double x0[N], double x1[N]);    void SORiteration2(double A[N][N], double b[N], double x0[N], double  x1[N]);    void triangularEquationsSolver(double T[N][N], double b[N], double x[N]);    int i, j;    static double A[N][N] = { { 5, 1, 1, 1, 1 }, { 1, 5, 1, 1, 1 }, { 1, 1, 5, 1, 1 }, { 1, 1, 1, 5, 1 }, { 1, 1, 1, 1, 5 } };    double b[N] = { 1.0, 1.0, 1.0, 1.0, 1.0 };    double x0[N] = { 1.0, 1.0, 1.0, 1.0, 1.0 };    double x00[N] = { 1.0, 1.0, 1.0, 1.0, 1.0 };    double xHf[N] = { 0 };    double x1[N];    do    {        memcpy(x0, x1, N*sizeof(double));        SORiteration1(A, b, x0, xHf);        SORiteration2(A, b, xHf, x1);    } while (normD(x0, x1) > epsilon);    printf("\n要求解的示例方程组为：\n A ||| b ||| x0\n");    for (i = 0; i < N; i++)    {        for (j = 0; j < N; j++)        {            printf("%f ", A[i][j]);        }        printf("||| %f||| %f\n", b[i], x00[i]);    }    printf("\n方程组解为：\n");    for (i = 0; i < N; i++)    {        printf("%f\n", x0[i]);    }    getchar();}double normD(double x0[N], double x1[N]){    int i;    double s = 0;    for (i = 0; i < N; i++)    {        s = s + (x0[i] - x1[i])*(x0[i] - x1[i]);    }    return sqrt(s);}void SORiteration1(double A[N][N], double b[N], double x0[N], double x1[N])//传数组往往是传其地址{    void triangularEquationsSolver(double T[N][N], double b[N], double x[N]);    //浪费一点内存没事，程序更具可读性。    double Lwl[N][N] = { 0 };    double Uwr[N][N] = { 0 };    //double bb[N];    int i, j;    double wb[N];    for (i = 1; i < N; i++)    {        for (j = 0; j < i; j++)        {            Lwl[i][j] = w*A[i][j];            Uwr[j][i] = -w*A[j][i];        }    }    for (i = 0; i < N; i++)    {        Lwl[i][i] = A[i][i];        Uwr[i][i] = (1 - w)*A[i][i];        wb[i] = w*b[i];    }    for (i = 0; i < N; i++)    {        for (j = i; j < N; j++)        {            wb[i] = wb[i] + Uwr[i][j] * x0[j];        }    }    triangularEquationsSolver(Lwl, wb, x1);    //printf("X1:\n");    //for (i = 0; i < N; i++)    //{    //  printf("%f", x1[i]);    //}    //getchar();}void SORiteration2(double A[N][N], double b[N], double x0[N], double x1[N]){    void triangularEquationsSolver(double T[N][N], double b[N], double x[N]);    //浪费一点内存没事，程序更具可读性。    double Uwl[N][N] = { 0 };    double Lwr[N][N] = { 0 };    int i, j;    double wb[N];    for (i = 1; i < N; i++)    {        for (j = 0; j < i; j++)        {            Uwl[j][i] = w*A[j][i];            Lwr[i][j] = -w*A[i][j];        }    }    for (i = 0; i < N; i++)    {        Uwl[i][i] = A[i][i];        Lwr[i][i] = (1 - w)*A[i][i];        wb[i] = w*b[i];        for (j = 0; j < i; j++)        {            wb[i] = wb[i] + Lwr[i][j] * x0[j];        }    }    triangularEquationsSolver(Uwl, wb, x1);}//求解上下三角方程的求解器void triangularEquationsSolver(double T[N][N], double b[N], double x[N]){    int i, j;    //for (i = 0; i < N; i++)    //{    //  for (j = 0; j < N; j++)    //  {    //      printf("%f  ", T[i][j]);    //  }    //  printf("\n");    //}    //getchar();    if (T[0][N - 1] == 0)    {        for (i = 0; i < N; i++)        {            for (j = 0; j < i; j++)            {                b[i] = b[i] - T[i][j] * x[j];            }            x[i] = b[i] / T[i][i];        }    }    else    {        for (i = N - 1; i >= 0; i--)        {            for (j = i + 1; j < N; j++)            {                b[i] = b[i] - T[i][j] * x[j];            }            x[i] = b[i] / T[i][i];        }    }}

  
  {width=”15cm”}\

• MATLAB代码

  clcclearw = 1;epsilon = 0.001;A = ones(5);for i = 1:5    A(i,i) = 5;endb = ones(5,1);x0 = b;D = diag(diag(A));CL = -tril(A - D);CU = -triu( A - D);Lwl = D - w*CL;Uwr = w*CU + (1-w)*D;Uwl = D - w*CU;Lwr = w*CL + (1-w)*D;wb = w*b;while (1)    xhf = Lwl\(Uwr*x0+wb);    x1 = Uwl \ (Lwr*xhf + wb);%   x1 = Lwl\(Uwr*x0+wb); %  x1 = (D-CL)^-1*CU*x0+(D-CL)^-1*b;   if(norm(x0-x1,2)<epsilon),break;end   x0 = x1;endx0A\bwb

各种迭代简单的比较

Jacobi迭代、Gauss-Seide迭代和最佳因子SOR迭代的比较，这里直接贴出结果：