求解三对角系统方程算法实现

算法名称：求解三对角系统方程

 

算法描述：

       线性方程系统的特例之一是三对角形式的，也就是说，非零的元素仅出现在对角线及其前、后一列的位置上。

       对三对角系统，LU分解过程、前代或回代每次只需O（N）次操作。

一般来说，不必把整个N×N矩阵全部保存起来，只需保存非零的元素，存为三向量。形如，

             （1）

注：a₀和c_N-1没有定义，示例程序中见不引用它们。

 

补充说明：

       在triadag中没有选主元。正是这个原因，即使矩阵是非奇异的，tridag也可能失败，因为有些非奇异矩阵是可能遇到零主元的。如果

                                                                              （2）

（称为对角占优）则可以证明算发不会遇到零主元。

       另外，可以构造出特殊的例子，来说明由于算法没有选主元而导致数值不稳定。然而，实践中这种不稳定几乎从未遇到过，不像一般的矩阵问题中选主元是必不可少的。

       一旦遇到tridag算法失败的问题，则可换用更一般的方法，即求解带状对角系统的方法。

 示例程序：

package  com.nc4nr.chapter02.tridag;

public   class  TriDag  ... {

    /**//*
     *        | 1.0, 2.0, 0.0,  0.0,   0.0 |
     *        | 3.0, 4.0, 5.0,  0.0,   0.0 |
     *    A = | 0.0, 6.0, 7.0,  8.0,   0.0 |
     *        | 0.0, 0.0, 9.0, 10.0,     11.0 |
     *        | 0.0, 0.0, 0.0, 12.0,    13.0 | 
     */
    
    double[]     b = ...{ 1.0, 4.0, 7.0, 10.0, 13.0 },
                a = ...{ 3.0, 6.0, 9.0, 12.0 },
                c = ...{ 2.0, 5.0, 8.0, 11.0 };
    
    /**//*        
     *         | 14.0 |
     *         | 15.0 |
     *    R = | 16.0 |
     *        | 17.0 |
     *        | 18.0 |
     */
    
    double[]    r = ...{ 14.0, 15.0, 16.0, 17.0, 18.0 };
    
    int anrow = 5;
    
    private void tridag() ...{
        double bet;
        int n = anrow;
        double[] gam = new double[n],
                 u = new double[n];
        if (b[0] == 0.0) System.out.println("Error 1 in tridag.");
        u[0] = r[0]/(bet = b[0]);
        for(int j = 1; j < n; j++) ...{
            gam[j] = c[j-1]/bet;
            bet = b[j]-a[j-1]*gam[j];
            if (bet == 0.0) System.out.println("Error 2 in tridag");
            u[j] = (r[j]-a[j-1]*u[j-1])/bet;
        }
        for (int j = n-2; j >= 0; j--)
            u[j] -= gam[j+1]*u[j+1];
        for (int j = 0; j < n; j++)
            System.out.println(u[j]);
    }
    
    public TriDag() ...{
        tridag();
    }
    
    public static void main(String[] args) ...{
        new TriDag();
    }

}