迭代法求解线性方程组

来源：互联网发布：量化对冲软件编辑：程序博客网时间：2024/04/28 11:40

关于迭代法一般解释参见数值计算教材
求解线性方程组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;    }

0 0