数值分析与多项式研究-Crout分解法

Crout分解即追赶法是专门用于求解三对角方程组的，这类方程组经常出现在用差分方法或有限元法求解二阶常微分方程边值问题、热传导问题及三次样条函数插值问题中，矩阵A的各阶顺序主子式不为0时有解，三对角矩阵各阶顺序主子式都不为0的一个充分条件是|a[i]|>|c[i]|+|b[i]|,c[i]*b[i]!=0,

即对角占优

当三对角矩阵满足对角占优时，追赶法是数值稳定的算法，具有程序简单，存贮量小和计算量小等优点。

#include <stdio.h>
#define N 1000
double a[N];
double b[N];
double c[N];
double d[N];
double x[N];
double y[N];
int n;
inline double abs(double x)
{if(x<0)return (-x);return x;}
int main()
{
    int i;
    int flag;
    while(scanf("%d",&n),n)
    {
        flag=1;
        for(i=0;i<n-1;++i)
            scanf("%lf",&c[i]);
        for(i=0;i<n;++i)
            scanf("%lf",&a[i]);
        for(i=1;i<n;++i)
            scanf("%lf",&b[i]);
        for(i=0;i<n;++i)
            scanf("%lf",&d[i]);
        for(i=0;i<n;++i)
        {
            if(i==0 && !(abs(a[i])>abs(c[i]) && c[i]!=0))
            {
                flag=0;
                break;
            }
            else if(i==(n-1) && !(abs(a[i])>abs(b[i]) && b[i]!=0))
            {
                flag=0;
                break;
            }
            else if(!(abs(a[i])>=(abs(c[i])+abs(b[i]))) || (c[i]*b[i]>-(1e-8) && c[i]*b[i]<1e-8))
            {
                flag=0;
                break;
            }
        }
        if(flag)
        {
            for(i=0;i<n;++i)
            {
                if(i!=0)
                    a[i]=a[i]-b[i]*c[i-1];
                if(a[i]==0)
                {
                    flag=0;
                    break;
                }
                c[i]=c[i]/a[i];
            }
        }
        if(flag)
        {
            y[0]=d[0]/a[0];
            for(i=1;i<n;++i)
                y[i]=(d[i]-y[i-1]*b[i])/a[i];
            x[n-1]=y[n-1];
            for(i=n-2;i>=0;--i)
                x[i]=y[i]-x[i+1]*c[i];
            for(i=0;i<n;++i)
                printf("%lf ",x[i]);
        }
        else
            printf("no unique solution");
        printf("\n");
    }
    return 0;
}