토마스 알고리즘에 의한 열 방정식의 해답

Question

토마스 알고리즘을 사용하여 미분 열 방정식을 풀려고합니다.

물리적 문제 : 플러그가 있고, 왼쪽에는 온도가 0, 오른쪽 온도는 1입니다.

토마스 알고리즘의 경우 방정식의 값 금액 인 QVector과 int을 허용하는 함수를 작성했습니다. 내가 컴파일하고 나는이 무엇입니까를 실행할 때 :

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 

또 다시 루프에서

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

내 문제는 다음과 같습니다

이 내 코드입니다.
나는 왜 이것을 얻고 있는지 전혀 모른다.

실수로 내 Thomas-method-function 결과가 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과 같습니다.

내가 첫 번째 단계는 주 대각선을 부정하는 것이었다 토마스 방법에 대해 읽어 기사에

, 내가 그것을했을, 그러나 함수의 끝에서 나는 그렇게 역 일을 수행해야합니다

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