title: Solution of heat equation by Thomas algorithm
body: <p>I am trying to solve differential heat equation using Thomas algorithm.  </p>
  
  <p><strong>Physical problem</strong>: We have plug, left side is having temperature <code>0</code>, right side temperature is <code>1</code>. </p>
  
  <p>For Thomas algorithm I have written a function, which accept three <code>QVector</code> and <code>int</code> value amount of equations.</p>
  
  <p>This is my code:</p>
  
  <pre><code>#include &lt;QCoreApplication&gt;
  #include &lt;QVector&gt;
  #include &lt;QDebug&gt;
  #include &lt;iostream&gt;
  using std::cin;
  
  void enterIn(QVector&lt;float&gt; &amp;Array, int Amount_of_elements){
      int transit;
      for(int i=0;i&lt;Amount_of_elements;i++){
          cin&gt;&gt;transit;
          Array.push_back(transit);
      }
  }
  
  QVector&lt;float&gt; shuttle_method(const QVector&lt;float&gt; &amp;below_main_diagonal,
                                QVector&lt;float&gt; &amp;main_diagonal,
                                const QVector&lt;float&gt; &amp;beyond_main_diagonal,
                                const QVector&lt;float&gt; &amp;free_term,
                                const int N){
      QVector &lt;float&gt; c;
      QVector &lt;float&gt; d;
  
      for(int i=0;i&lt;N;i++){
          main_diagonal[i]*=(-1);
      }
  
      QVector&lt;float&gt; 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&lt;=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()&lt;&lt;x.size()&lt;&lt;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&gt;=0;i--){
          x[i]=c[i]*x[i+1]+d[i];
          // qDebug()&lt;&lt;x[i]&lt;&lt;endl;
      }
  
      return x;
  }
  
  int main()
  {    
      QVector &lt;float&gt; alpha;  // below
      QVector &lt;float&gt; beta;   // main diagonal * (-1)
      QVector &lt;float&gt; gamma;  // beyond
      QVector &lt;float&gt; b;      // free term
  
      QVector&lt;float&gt; 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&lt;cells_x;i++){
          alpha[i]= -dt/(h*h);
      }
  
      for(int i=0;i&lt;cells_x;i++){
          beta[i] = (2*dt)/(h*h)+1;
      }
      for(int i=0;i&lt;cells_x-1;i++){
          gamma[i]= -dt/(h*h);
      }
      gamma[cells_x-1]=0;
  
      qDebug()&lt;&lt;"alpha= "&lt;&lt;endl&lt;&lt;alpha.size()&lt;&lt;alpha&lt;&lt;endl&lt;&lt;"beta = "&lt;&lt;endl&lt;&lt;beta.size()&lt;&lt;beta&lt;&lt;endl&lt;&lt;"gamma= "&lt;&lt;gamma.size()&lt;&lt;gamma&lt;&lt;endl;
  
  
      for(int i=0;i&lt;cells_x-1;i++){
          T[i]=0;
      }
      T[cells_x-1]=1;
      qDebug()&lt;&lt;endl&lt;&lt;endl&lt;&lt;T&lt;&lt;endl;
  
      //qDebug()&lt;&lt; shuttle_method(alpha,beta,gamma,b,N);
  
      QVector&lt;float&gt; 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 &lt; 50; stepTime++){
          Tn = shuttle_method(alpha,beta,gamma,Tn,cells_x);
          Tn[0]=0;
          Tn[cells_x-1]=1;
          qDebug()&lt;&lt;Tn&lt;&lt;endl;
      }
      return 0;
  }
  </code></pre>
  
  <p>My problem is:
  when I compile and run it I am getting this:</p>
  
  <pre><code>Tn  &lt;20 items&gt;  QVector&lt;float&gt;
  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  &lt;20 items&gt;  QVector&lt;float&gt;
  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  &lt;20 items&gt;  QVector&lt;float&gt;
  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  &lt;20 items&gt;  QVector&lt;float&gt;
  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
  </code></pre>
  
  <p>Again and again in loop.<br>
  I have no idea why I am getting this. </p>
  
  <p>Maybe my mistake is in assign <code>Tn</code> result of my Thomas-method-function?<br>
  or in realization of Thomas method? or in boundary conditions?</p>