线性方程组解的迭代改进算法实现

算法名称：

线性方程组解的迭代改进

 

算法描述：

       假定向量x是线性方程组

                                                                             （1）

的精确解，然而并不知道解x，而仅仅知道一些略有误差的解x+δx，其中δx为未知的误差。当将其乘以矩阵A时，这个有误差的解使乘积与原来的右端项b也存在偏差，即

                                                                            （2）

由（2）式减去（1）式得，

                                                                                          （3）

但由（2）式也可求得δb，将其带入(3)式中得，

                               （4）

该方程中整个右端项是已知量，因为x+δx是有待改进的有误差的解。然后，只需求解（4）式的偏差量δx，将其从有误差的解中减去便得到改进后的解。

       这种情况下，已经有了矩阵A的LU分解形式，要求解（4）式，只需计算右端项再进行回代即可。

运行示例：

The origin matrix:

0.0 3.0 1.0 1.0

2.0 2.0 3.0 2.0

4.0 -3.0 0.0 1.0

6.0 1.0 -2.0 -5.0

--------------------------

The matrix (after LU decomposition):

6.0 1.0 -2.0 -5.0

0.0 3.0 1.0 1.0

0.3333333333333333 0.5555555555555556 3.1111111111111107 3.1111111111111107

0.6666666666666666 -1.222222222222222 0.8214285714285714 3.0000000000000004

--------------------------

The right vector:

0.0

-2.0

-5.0

6.0

--------------------------

The solution vector (before improvement):

-0.4523809523809522

0.42857142857142866

0.6190476190476184

-1.9047619047619044

--------------------------

The solution vector (after improvement):

-0.45238095238095233

0.4285714285714286

0.6190476190476188

-1.9047619047619049

--------------------------

 

代码示例：

package  com.nc4nr.chapter02.mprove;

public class Mprove ... {

    // 建议从构造函数看起，Mprove
    
    // 4 * 4 coefficient matrix a
    double[][] a = ...{
            ...{0.0, 3.0, 1.0, 1.0},
            ...{2.0, 2.0, 3.0, 2.0},
            ...{4.0, -3.0, 0.0, 1.0},
            ...{6.0, 1.0, -2.0, -5.0}
    };
    
    // 4 * 1 coefficient matrix b
    double[][] b = ...{
            ...{0.0},
            ...{-2.0},
            ...{-5.0},
            ...{6.0}
    };
    
    int anrow = 4;
    int[] indx = new int[anrow];
    double[][] x = new double[anrow][1]; // 请留意后面的注释
    double[][] ora = new double[anrow][anrow];    // 拷贝分解前的系数矩阵a，因为a被分解后将被破坏
    
    private void ludcmp() ...{
        int n = anrow, imax = 0;
        double[] vv = new double[anrow];
        for (int i = 0; i < n; i++) ...{
            double big = 0.0;
            for (int j = 0; j < n; j++) ...{
                double temp = Math.abs(a[i][j]);
                if (temp > big) big = temp;
            }
            if (big == 0.0) System.out.println("lu: singular matrix");
            vv[i] = 1.0/big;
        }
        
        for (int j = 0; j < n; j++) ...{
            for (int i = 0; i < j; i++) ...{
                double sum = a[i][j];
                for (int k = 0; k < i; k++) sum -= a[i][k]*a[k][j];
                a[i][j] = sum;
            }
            double big = 0.0;
            for (int i = j; i < n; i++) ...{
                double sum = a[i][j];
                for (int k = 0; k < j; k++) sum -= a[i][k]*a[k][j];
                a[i][j] = sum;
                double dum = Math.abs(vv[i]*sum);
                if (dum >= big) ...{
                    big = dum;
                    imax = i;
                }
            }
            if (j != imax) ...{
                for (int i = 0; i < n; i++) ...{
                    double mid = a[imax][i];
                    a[imax][i] = a[j][i];
                    a[j][i] = mid;
                }
                vv[imax] = vv[j];
            }
            indx[j] = imax;
            if (j != n-1) ...{
                double dum = 1.0/a[j][j];
                for (int k = j+1; k < n; k++) a[k][j] *= dum;
            }
        }
    }
    
    private void lubksb(double[][] alu, int[] idx, double[][] rb) ...{
        int n = anrow, ii = 0;
        for (int i = 0; i < n; i++) ...{
            int k = idx[i];
            double sum = rb[k][0];
            rb[k][0] = rb[i][0];
            if (ii != 0)
                for (int j = 0; j < i; j++) sum -= alu[i][j]*rb[j][0];
            else if (sum != 0)
                ii = i+1;
            rb[i][0] = sum;
        }
        for (int i = n-1; i >= 0; i--) ...{
            double dum = 1.0/alu[i][i];
            double sum = rb[i][0];
            for (int k = i+1; k < n; k++) sum -= alu[i][k]*rb[k][0];
            rb[i][0] = dum * sum;
        }
    }
    
    private void mprove() ...{
        int n = anrow;
        double[][] r = new double[anrow][1];
        for (int i = 0; i < n; i++) ...{
            double sdp = -b[i][0];
            for (int j = 0; j < n; j++) ...{
                sdp += ora[i][j]*x[j][0];
            }
            r[i][0] = sdp;
        }
        lubksb(a,indx,r);
        for (int i = 0; i < n; i++) x[i][0] -= r[i][0];
    }
    
    private void output(double[][] m, int nrow, int ncol) ...{
        for (int i = 0; i < nrow; i++) ...{
            String str = "";
            for (int j = 0; j < ncol; j++) str += m[i][j] + " ";
            System.out.println(str);
        }
        System.out.println("--------------------------");
    }
    
    public Mprove() ...{
        // 拷贝原始系数矩阵a到ora
        for (int i = 0; i < anrow; i++) ...{
            for (int j = 0; j < anrow; j++) ...{
                ora[i][j] = a[i][j];
            }
        }
        // 拷贝原始右端项b
        for (int i = 0; i < anrow; i++) ...{
            x[i][0] = b[i][0];
        }
        System.out.println("The origin matrix:");
        output(ora,anrow,anrow);
        // 这里将对a进行LU分解，并将分解结果存放到a
        ludcmp();
        System.out.println("The matrix (after LU decomposition):");
        output(a,anrow,anrow);
        System.out.println("The right vector:");
        output(b,anrow,1);
        // 这里传入的是LU分解后的矩阵a，x将为带误差的解向量
        lubksb(a,indx,x);
        System.out.println("The solution vector (before improvement):");
        output(x,anrow,1);
        // 这里将进行迭代改进x
        mprove();
        System.out.println("The solution vector (after improvement):");
        output(x,anrow,1);
        // 本程序没有改变b和ora
    }
    
    public static void main(String[] args) ...{
        new Mprove();
    }

}