迭代法求解线性方程组

关于迭代法一般解释参见数值计算教材
求解线性方程组 Ax=b

Jacobi 迭代法

向量 Xk+1=D−1[−(L+U)Xk+b]
其中 D=diag[a11,a12,...,ann]T
L 为系数矩阵的下三角矩阵，其中除掉对角线
U 为系数矩阵的上三角矩阵，除去对角线

java代码

/** *Jacobi迭代法Ax = b,若在给定步数能求出解，返回true，否则返回false * @param A * @param b * @param x解 * @param n未知数个数 * @param esp误差 */
    public static   boolean  Jacobi(double[][] A,double[] b,double[]x,int n,double esp,int step)
    {
        double[] y = new double[n+1];
        for(int i=1 ; i<=n ; ++i)x[i] = 0;
        int k = 0;
        boolean flag = false;
        while(k<step)
        {
            for(int i=1 ; i<=n ; ++i)
            {
                double sum = b[i];
                for(int j=1 ; j<=n ; ++j)
                {
                    if(i!=j)
                    {
                        sum-=A[i][j]*x[j];
                    }
                }
                y[i] = sum/A[i][i];
            }

            //比较误差
            double max = Math.abs(y[0]-x[0]);
            for(int i=2 ; i<=n ; ++i)
                if(max>esp)
                {
                    break;
                }else{
                    max = Math.max(max, Math.abs(y[i]-x[i]));
                }
            for(int i=1 ; i<=n ; ++i)x[i] = y[i];
            if(max<esp){
                flag = true;
                break;
            }
            k++;
        }
        return flag;
    }

空间复杂度 O(n2) ,时间复杂度 O(k∗n2) ,该算法 还可以优化

Gauss_Seidel

这种方法只是在 i<j 的时候用于计算新向量的值用最新算出的解来代替，大同小异

    /** * 求解Ax =b * @param A * @param b * @param x * @param step步数 * @param n未知元个数 * @param exp误差限 * @return true IFF 在给定步数内求出解 */
    public static boolean GauusSeidel(double[][] A,double[]b,double[]x,int step,int n,double exp)
    {
        double[] y = new double[n+1];
        for(int i=1 ; i<=n ; ++i)x[i] = 0;
        int k = 0;
        boolean flag = false;
        while(k<step)
        {
            for(int i=1 ; i<=n ; ++i)
            {
                double sum = b[i];
                for(int j=1 ; j<=n ; ++j)
                {
                    if(i>j)
                    {
                        sum-=A[i][j]*y[j];
                    }else if(i<j)
                    {
                        sum-=A[i][j]*x[j];
                    }
                }
                y[i] = sum/A[i][i];
            }

            //比较误差
            double max = Math.abs(y[0]-x[0]);
            for(int i=2 ; i<=n ; ++i)
                if(max>exp)
                {
                    break;
                }else{
                    max = Math.max(max, Math.abs(y[i]-x[i]));
                }
            for(int i=1 ; i<=n ; ++i)x[i] = y[i];
            if(max<exp){
                flag = true;
                break;
            }
            k++;
        }
        return flag;
    }