数值分析，高斯赛德尔，最速下降，共轭梯度迭代法求解

一、题目要求

二、原理分析

2.1 高斯-赛德尔方法原理

设方程组Ax = f 有唯一解，将Ax = b 变形为等价的方程组 x = Bx + b，将公式化简为以下形式：

                   j=1naijxj= bi(i=1,2…n)                  （式1）

由此可以建立迭代公式，如式二所示：

 xi(k+1)=1aii[bi-j=1i-1aijxj(k+1)-j=i+1naijxj(k)] (i=1,2…n;k=1,2,……)    （式2）

取初试向量x(0)=[x10,x20…xn(0)]T，进行迭代计算，高斯-赛德尔迭代法比简单迭代法收敛快，相比简单迭代法，高斯-赛德尔跌打法在计算xik+1时利用了刚刚迭代出的x1k+1，x2k+1，…，xi-1k+1。

              收敛性：已知定理，如果系数矩阵A对称正定，则高斯赛德尔迭代公式收敛。已知系数矩阵A是Hilbert矩阵，是对称正定的，所以高斯赛德尔迭代公式收敛。

2.2 最速下降法原理

       已知A =（aii）n*n为实对称正定矩阵，b，x∈Rn，则线性方程组Ax = b的解x满足：解x使得二次函数φx= 12Ax,x-(b,x)取得极小值。所以求解线性方程组Ax = b的解，即为求解二次函数φx= 12Ax,x-(b,x)的极小值点。求二次函数φx极小值点的一般方法是：构造一个向量序列{x(k)},使得φxk→minφx。

方法步骤如下：

       1）任取一个初始向量x0；

       2）构造迭代格式xk+1=xk+αkpk，其中pk是搜索方向，αk是搜索步长；

       3）选择pk和αk使得φxk+1 = φ（xk+αkpk）<φxk,则当k→∞时，有φxk→φx* = minφx；

       4）计算误差ε，直到||xk-xk-1||<ε或者||rk||<ε迭代结束。

       在以上步骤中，对于搜索步长αk，当αk = rk，pkApk，pk时，αk是k到k+1的最优步长。对于搜索方向pk取为函数φx减少最快的方向，即φx的负梯度方向，此时pk=rk。

     最速下降法就是使得目标函数下降的最快的方向，反映了目标函数的一种局部性质，最速下降法是线性收敛的算法。

       收敛性：最优步长规则下的最速下降法是线性收敛的算法。

2.3 共轭梯度法原理

       共轭梯度法同时考量上一步和当前方向，目前的下降方向是上步下降方向和当前残量的线性组合。pk=rk+βk-1pk-1。寻求步长α，使得φx+αp = min，其中αk= rk，pkApk，pk，对于pk=rk+βk-1pk-1，选择βk-1使得φx+αp = min，可得βk-1= pk-1TArkpk-1TApk-1，综合以上概念可得到共轭梯度发的基本流程。

       1）任取一个初始向量x（0） ；

       2）构造迭代格式xk+1=xk+αkpk，其中pk是搜索方向，αk是搜索步长；

       3）计算步长αk = rk，rkApk，pk；

       4）rk+1 = rk - αkApk；

       5）βk= rk+1，rk+1rk，rk；

       6）pk= rk+1+βkpk；

       7）计算误差ε，直到||rk||<ε迭代结束。

收敛性：最优步长规则下的共轭梯度法是线性收敛的算法。

三、代码实现

 代码用java语言实现

public class Main {
    public static void main(String[] args) {
/**
 * 数值分析作业
 */

        double[][] A = hilbert(12);//设A为一个12阶次的hilbert矩阵
        double[] b = {1353.0/436,1295.0/594,1121.0/640,934.0/629,2142.0/1651,1162.0/1005
        ,1645.0/1574,1545.0/1618,901.0/1024,787.0/964,611.0/802,793.0/1110};
        double[] x0 = new double[12];
        GS(A,b,x0);//高斯赛德尔方法
        GD(A,b,x0);//梯度下降法（x-x0）
        GDR(A,b,x0);//梯度下降法（r）
        CG(A,b,x0);//共轭梯度法


    }

    static double[][] hilbert(int n){//定义hilbert矩阵的方法
        double[][] A = new double[n][n];
        for (int i = 0; i < A.length; i++) {
            for (int j = 0; j < A[i].length; j++) {
                A[i][j] = (1.0/(i+j+1));
            }
        }
        return A;
    }


    static void GS(double[][] A,double[] b,double[] x0){//高斯-赛德尔迭代法
        double[] x =Arrays.copyOf(x0,x0.length);
        int number = 0;

        while(true){
            for (int i = 0; i < A.length; i++) {//计算迭代后的下一次x
                double temp = 0;
                for (int j = 0; j < A[i].length; j++) {
                    if (i != j) {
                        temp += x[j] * A[i][j];
                    }
                }
                x[i] = (b[i] - temp) / A[i][i];
            }
            number++;
            double[] err = new double[x.length];
            for(int i = 0; i< x.length; i++){//以x-x0的无穷范数为判断标准
                err[i] = Math.abs(x[i] - x0[i]);
            }
            Arrays.sort(err);
            if(err[err.length-1] < 0.0001){//当无穷范数小于误差值输出结果
                System.out.println(String.format("高斯-赛德尔迭代法：迭代次数k =%d,其X值为：",number)
                        +Arrays.toString(x));
               break;
            }
            x0 = Arrays.copyOf(x,x.length);
        }
    }

    static void GD(double[][] A,double[] b,double[] x0){//最速下降法方法
        double[] x =Arrays.copyOf(x0,x0.length);
        int number = 0;
        double[] r = new double[x0.length];
        while(true){
            for(int i = 0; i<A.length;i++){//计算迭代后的r
                double temp = 0;
                for(int j = 0 ;j < A[i].length;j++){
                    temp += x[j] * A[i][j];
                }
                r[i] = b[i] - temp;
            }
            double element = 0 ;
            for(int i = 0;i<r.length;i++){//计算步长式子的分子项
                element += r[i]*r[i];
            }
            double[] r1 = new double[r.length];
            for(int i = 0; i<A.length;i++){
                double temp = 0;
                for(int j = 0 ;j < A[i].length;j++){
                    temp += r[j] * A[i][j];
                }
                r1[i] = temp;
            }
            double denominator = 0;
            for(int i = 0;i<r.length;i++) {//计算步长式子的分母项
                denominator += r1[i] * r[i];
            }
            double a = element/denominator;//计算步长
            for(int i = 0;i<r.length;i++){
                x[i] = x[i]+a*r[i];
            }
            number++;

            double[] err = new double[x.length];
            for(int i = 0; i< x.length; i++){//计算x-x0的绝对值为err值，即无穷范数
                err[i] = Math.abs(x[i] - x0[i]);
            }
            Arrays.sort(err);
            if(err[err.length-1] < 0.0001){//当x-x0的绝对范数小于误差
                System.out.println(String.format("最速下降迭代法（x-x0）：迭代次数k =%d,其X值为：",number)
                        +Arrays.toString(x));
                break;
            }
            x0 = Arrays.copyOf(x,x.length);

        }
    }

    static void GDR(double[][] A,double[] b,double[] x0){//最速下降法方法
        double[] x =Arrays.copyOf(x0,x0.length);
        int number = 0;
        double[] r = new double[x0.length];
        while(true){
            for(int i = 0; i<A.length;i++){//计算迭代后的r
                double temp = 0;
                for(int j = 0 ;j < A[i].length;j++){
                    temp += x[j] * A[i][j];
                }
                r[i] = b[i] - temp;
            }
            double element = 0 ;
            for(int i = 0;i<r.length;i++){//计算步长式子的分子项
                element += r[i]*r[i];
            }
            double[] r1 = new double[r.length];
            for(int i = 0; i<A.length;i++){
                double temp = 0;
                for(int j = 0 ;j < A[i].length;j++){
                    temp += r[j] * A[i][j];
                }
                r1[i] = temp;
            }
            double denominator = 0;
            for(int i = 0;i<r.length;i++) {//计算步长式子的分母项
                denominator += r1[i] * r[i];
            }
            double a = element/denominator;//计算步长
            for(int i = 0;i<r.length;i++){
                x[i] = x[i]+a*r[i];
            }
            number++;

            double[] err = new double[x.length];
            for(int i = 0; i< x.length; i++){//计算r的绝对值为err值，即无穷范数
                err[i] = Math.abs(r[i]);
            }
            Arrays.sort(err);
            if(err[err.length-1] < 0.0001){//当r的绝对范数小于误差
                System.out.println(String.format("最速下降迭代法（r）：迭代次数k =%d,其X值为：",number)
                        +Arrays.toString(x));
                break;
            }
            x0 = Arrays.copyOf(x,x.length);

        }
    }

    static void CG(double[][] A,double[] b,double[] x0){//共轭梯度法方法
        double[] x =Arrays.copyOf(x0,x0.length);//复制内容
        int number = 0;
        double[] r = new double[x0.length];
        for(int i = 0; i<A.length;i++){//求r
            double temp = 0;
            for(int j = 0 ;j < A[i].length;j++){
                temp += x[j] * A[i][j];
            }
            r[i] = b[i] - temp;
        }
        double[] p = Arrays.copyOf(r,r.length);//p(0) = r(0) = b
        double[] ap = new double[r.length];
        while(true){
            double[] r0 = Arrays.copyOf(r,r.length);//保留本次r值
            for(int i = 0; i<A.length;i++){//求AP
                double temp = 0;
                for(int j = 0 ;j < A[i].length;j++){
                    temp += p[j] * A[i][j];
                }
                ap[i] = temp;
            }
            double element = 0 ;
            for(int i = 0;i<r.length;i++){//求r与r的内积
                element += r[i]*r[i];
            }
            double denominator = 0;
            for(int i = 0;i<r.length;i++) {//求ap与p的内积
                denominator += p[i] * ap[i];
            }
            double a = element/denominator;//求a
            for(int i = 0;i<r.length;i++){//求x迭代一次的值
                x[i] = x[i]+a*p[i];
            }
            number++;//记录迭代次数
            for(int i = 0; i<A.length;i++){//求下一次的r
                r[i] = r[i] - a*ap[i];
            }
            double[] err = new double[x.length];
            for(int i = 0; i< x.length; i++){//计算梯度r的无穷范数即无穷范数
                err[i] = Math.abs(r[i]);
            }
            Arrays.sort(err);
            if(err[err.length-1]<0.0001){//终止迭代条件，梯度趋于0
                System.out.println(String.format("共轭梯度迭代法：迭代次数k =%d,其X值为：",number)
                        +Arrays.toString(x));
                break;
            }
            x0 = Arrays.copyOf(x,x.length);


            double element1 = 0 ;
            for(int i = 0;i<r.length;i++){//求r(k+1)与r(k+1)的内积
                element1 += r[i]*r[i];
            }
            double denominator1 = 0;
            for(int i = 0;i<r.length;i++) {//求r(k)与r(k)的内积
                denominator1 += r0[i] * r0[i];
            }
            double beta = element1/denominator1;//求beta
            for(int i = 0; i<A.length;i++){//求下一次的p
                p[i] = r[i] + beta*p[i];
            }
        }
    }

}

四、结果展示

最速下降法用||xk-xk-1||<ε和||rk||<ε两种方式判断迭代结束。