Thomas算法求解熱方程

我想用Thomas算法求解差熱方程。Thomas算法求解熱方程

身體問題：我們有插頭，左邊有溫度0，右邊的溫度是1。

對於托馬斯算法我寫了一個函數，它接受三個方程的值爲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

連連。
我不知道爲什麼我得到這個。

也許我的錯誤是分配Tn我的托馬斯方法函數的結果？
還是在實現托馬斯方法？或在邊界條件？

來源

2016-09-27 John

您是否嘗試過運行和步進通過調試器中的代碼？這通常有幫助。從小號開始，這樣循環就能快速完成。 –

解決這些問題的正確工具是您的調試器。在*堆棧溢出問題之前，您應該逐行執行您的代碼。如需更多幫助，請閱讀[如何調試小程序（由Eric Lippert撰寫）]（https://ericlippert.com/2014/03/05/how-to-debug-small-programs/）。至少，您應該\編輯您的問題，以包含一個[最小，完整和可驗證]（http://stackoverflow.com/help/mcve）示例，該示例再現了您的問題，以及您在調試器。 –

當然，我用它！但它沒有幫助，我仍然得到相同的結果，並不能得到爲什麼 – John

我明白了！

邊界條件必須採取行動矢量

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]將總是等於零並且T [N-1]將等於1;

而且在那裏我讀到托馬斯方法的第一步是要否定主對角線上的文章中，我已經做到了，但隨後在函數結束時，我必須做相反的事情，所以：

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; 
    }

，並在最後一步，我們會得到絕對的「物理」的結果：

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

來源

2016-09-29 15:06:12 John

Thomas算法求解熱方程

回答

相關問題