Thomas算法求解热方程

Question

我正在尝试使用 Thomas 算法求解微分热方程。

物理问题：我们有插头，左边温度0，右边温度1。

对于 Thomas 算法，我写了一个函数，它接受三个 QVector 和 int 值的方程。

这是我的代码：

#include <QCoreApplication>
#include <QVector>
#include <QDebug>
#include <iostream>
using std::cin;

void enterIn(QVector<float> &Array, int Amount_of_elements){
    int transit;
    for(int i=0;i<Amount_of_elements;i++){
        cin>>transit;
        Array.push_back(transit);
    }
}

QVector<float> shuttle_method(const QVector<float> &below_main_diagonal,
                              QVector<float> &main_diagonal,
                              const QVector<float> &beyond_main_diagonal,
                              const QVector<float> &free_term,
                              const int N){
    QVector <float> c;
    QVector <float> d;

    for(int i=0;i<N;i++){
        main_diagonal[i]*=(-1);
    }

    QVector<float> x; //result

    c.push_back(beyond_main_diagonal[0]/main_diagonal[0]);
    d.push_back(-free_term[0]/main_diagonal[0]);

    for(int i=1;i<=N-2;i++){
        c.push_back(beyond_main_diagonal[i]/(main_diagonal[i]-below_main_diagonal[i]*c[i-1]));
        d.push_back(  (below_main_diagonal[i]*d[i-1]  -  free_term[i])  /  (main_diagonal[i]-  below_main_diagonal[i]*c[i-1])  );
    }

    x.resize(N);
    //qDebug()<<x.size()<<endl;

    int n=N-1;
    x[n]=(below_main_diagonal[n]*d[n-1]-free_term[n])/(main_diagonal[n]-below_main_diagonal[n]*c[n-1]);

    for(int i=n-1;i>=0;i--){
        x[i]=c[i]*x[i+1]+d[i];
        // qDebug()<<x[i]<<endl;
    }

    return x;
}

int main()
{    
    QVector <float> alpha;  // below
    QVector <float> beta;   // main diagonal * (-1)
    QVector <float> gamma;  // beyond
    QVector <float> b;      // free term

    QVector<float> T;

    int cells_x=40;         //amount of equations
    alpha.resize(cells_x);
    beta.resize(cells_x);
    gamma.resize(cells_x);
    b.resize(cells_x);
    T.resize(cells_x);

    float dt=0.2,h=0.1;

    alpha[0]=0;
    for(int i=1;i<cells_x;i++){
        alpha[i]= -dt/(h*h);
    }

    for(int i=0;i<cells_x;i++){
        beta[i] = (2*dt)/(h*h)+1;
    }
    for(int i=0;i<cells_x-1;i++){
        gamma[i]= -dt/(h*h);
    }
    gamma[cells_x-1]=0;

    qDebug()<<"alpha= "<<endl<<alpha.size()<<alpha<<endl<<"beta = "<<endl<<beta.size()<<beta<<endl<<"gamma= "<<gamma.size()<<gamma<<endl;


    for(int i=0;i<cells_x-1;i++){
        T[i]=0;
    }
    T[cells_x-1]=1;
    qDebug()<<endl<<endl<<T<<endl;

    //qDebug()<< shuttle_method(alpha,beta,gamma,b,N);

    QVector<float> Tn;
    Tn.resize(cells_x);
    Tn = shuttle_method(alpha,beta,gamma,T,cells_x);
    Tn[0]=0;Tn[cells_x-1]=1;
    for(int stepTime = 0; stepTime < 50; stepTime++){
        Tn = shuttle_method(alpha,beta,gamma,Tn,cells_x);
        Tn[0]=0;
        Tn[cells_x-1]=1;
        qDebug()<<Tn<<endl;
    }
    return 0;
}

我的问题是：当我编译并运行它时，我得到这个：

Tn  <20 items>  QVector<float>
0   float
0.000425464 float
0.000664658 float
0.000937085 float
0.00125637  float
0.00163846  float
0.00210249  float
0.00267163  float
0.00337436  float
0.00424581  float
0.00532955  float
0.00667976  float
0.00836396  float
0.0104664   float
0.0130921   float
0.0163724   float
0.0204714   float
0.0255939   float
0.0319961   float

Tn  <20 items>  QVector<float>
0   float
-0.000425464    float
0.000643385 float
-0.000926707    float
0.00120951  float
-0.00161561 float
0.00202056  float
-0.00263167 float
0.00324078  float
-0.00418065 float
0.00511726  float
-0.00657621 float
0.00802998  float
-0.0103034  float
0.0125688   float
-0.0161171  float
0.0196527   float
-0.0251945  float
0.0307164   float
1   float

Tn  <20 items>  QVector<float>
0   float
0.000425464 float
0.000664658 float
0.000937085 float
0.00125637  float
0.00163846  float
0.00210249  float
0.00267163  float
0.00337436  float
0.00424581  float
0.00532955  float
0.00667976  float
0.00836396  float
0.0104664   float
0.0130921   float
0.0163724   float
0.0204714   float
0.0255939   float
0.0319961   float

Tn  <20 items>  QVector<float>
0   float
-0.000425464    float
0.000643385 float
-0.000926707    float
0.00120951  float
-0.00161561 float
0.00202056  float
-0.00263167 float
0.00324078  float
-0.00418065 float
0.00511726  float
-0.00657621 float
0.00802998  float
-0.0103034  float
0.0125688   float
-0.0161171  float
0.0196527   float
-0.0251945  float
0.0307164   float
1   float

周而复始。
我不知道为什么我会得到这个。

也许我的错误在于我的 Thomas 方法函数的赋值 Tn 结果？
或实现托马斯方法？还是在边界条件下？

Answer 1

我明白了！

边界条件必须作用于向量

QVector<float> below_main_diagonal, 
QVector<float>  main_diagonal, 
QVector<float>  beyond_main_diagonal

所以T[0]一定是0，T[N-1]一定是1。我们可以这样：

main_diagonal.first()=1;
main_diagonal.last()=1;
beyond_main_diagonal.first()=0;
below_main_diagonal.last()=0;

因此，T[0] 将始终等于 0，而 T[N-1] 将始终等于 1；

在我读到 Thomas 方法的文章中，第一步是取反主对角线，我已经做到了，但是在函数的最后我必须做相反的事情，所以：

for(int i(0);i<N;++i){
   main_diagonal[i]*=(-1);
}

我们可以再次使用这个功能，这不是绝对最优，但它工作稳定。

然后，整个代码将如下所示：

#include <QCoreApplication>
#include <QVector>
#include <QDebug>
#include <iostream>


QVector<float> Thomas_Algorithm( QVector<float> &below_main_diagonal ,
                                     QVector<float> &main_diagonal  ,
                                     QVector<float> &beyond_main_diagonal ,
                                     QVector<float> &free_term,
                                     const  int N){

        QVector<float> x; //vector of result

        // checking of input data
        if(below_main_diagonal.size()!=main_diagonal.size()||
                main_diagonal.size()!=beyond_main_diagonal.size()||
                free_term.size()!=main_diagonal.size())
        { qDebug()<<"Error!\n"
                    "Error with accepting Arrays! Dimensities are different!"<<endl;
            x.resize(0);
            return x;
        }
        if(below_main_diagonal[0]!=0){
            qDebug()<< "Error!\n"
                       "First element of below_main_diagonal must be equal to zero!"<<endl;
            x.resize(0);
            return x;
        }
        if(beyond_main_diagonal.last()!=0){
            qDebug()<< "Error!\n"
                       "Last element of beyond_main_diagonal must be equal to zero!"<<endl;
            x.resize(0);
            return x;

        }
        // end of checking

        QVector <float> c;
        QVector <float> d;

        for(int i=0;i<N;i++){
            main_diagonal[i]*=(-1);
        }

        c.push_back(beyond_main_diagonal[0]/main_diagonal[0]);
        d.push_back(-free_term[0]/main_diagonal[0]);

        for(int i=1;i<=N-2;i++){
            c.push_back(beyond_main_diagonal[i]/(main_diagonal[i]-below_main_diagonal[i]*c[i-1]));
            d.push_back(  (below_main_diagonal[i]*d[i-1]  -  free_term[i])  /
                    (main_diagonal[i]-  below_main_diagonal[i]*c[i-1])  );
        }

        x.resize(N);

        int n=N-1;
        x[n]=(below_main_diagonal[n]*d[n-1]-free_term[n])/(main_diagonal[n]-below_main_diagonal[n]*c[n-1]);

        for(int i=n-1;i>=0;i--){
            x[i]=c[i]*x[i+1]+d[i];
        }

        for(int i(0);i<N;++i){
           main_diagonal[i]*=(-1);
        }

        return x;
    }

    int main()
    { 
        QVector <float> alpha;  // below
        QVector <float> beta;   // main diagonal * (-1)
        QVector <float> gamma;  // beyond
        QVector <float> b;      // free term


        QVector<float> T;

        int cells_x=30;         // amount of steps
        alpha.resize(cells_x);
        beta.resize(cells_x);
        gamma.resize(cells_x);
        T.resize(cells_x );

        float dt=0.2,h=0.1;

        alpha[0]=0;
        for(int i=1;i<cells_x-1;i++){
            alpha[i]= -dt/(h*h);
        }
        alpha[cells_x-1]=0;

        beta[0]=1;
        for(int i=1;i<cells_x-1;i++){
            beta[i] = (2*dt)/(h*h)+1;
        }
        beta[cells_x-1]=1;
        gamma[0]=0;
        for(int i=1;i<cells_x-1;i++){
            gamma[i]= -dt/(h*h);
        }
        gamma[cells_x-1]=0;  

        for(int i=0;i<cells_x-1;i++){
            T[i]=0;
        }
        T[cells_x-1]=1;

        QVector<float>Tn;
        Tn.resize(cells_x);
        Tn=  Thomas_Algorithm(alpha,beta,gamma,T,cells_x);
        // boundary conditions!
        beta.first()=1;
        beta.last()=1;
        gamma.first()=0;
        alpha.last()=0;
        // and then due to bc we always have T[0]=0 and T[n]=1

        for(int stepTime=0;stepTime<100;stepTime++){
            Tn = Thomas_Algorithm(alpha,beta,gamma,Tn,cells_x);

            qDebug()<<"stepTime = "<<stepTime<<endl<<endl;
            qDebug()<<Tn<<endl;

            // boundary conditions!
            beta.first()=1;
            beta.last()=1;
            gamma.first()=0;
            alpha.last()=0;
            // and then due to bc we always have T[0]=0 and T[n]=1
        }

        return 0;
    }

在最后一步中，我们绝对会得到 "physical" 结果：

Tn  <30 items>  QVector<float>
                0   float
                0.0344828   float
                0.0689656   float
                0.103448    float
                0.137931    float
                0.172414    float
                0.206897    float
                0.24138 float
                0.275862    float
                0.310345    float
                0.344828    float
                0.379311    float
                0.413793    float
                0.448276    float
                0.482759    float
                0.517242    float
                0.551724    float
                0.586207    float
                0.62069 float
                0.655173    float
                0.689655    float
                0.724138    float
                0.758621    float
                0.793104    float
                0.827586    float
                0.862069    float
                0.896552    float
                0.931035    float
                0.965517    float
                1   float

Thomas算法求解热方程

Solution of heat equation by Thomas algorithm

c++

equation

heat