この2つのコードの違いは何ですか？

Question

差分問題を解くための2つのコードを書きました。どちらも予測子補正ツールを使っています。

最初は差分系の解法でもありますので、基本的には系で解く方程式の数に等しい大きさのベクトルyを使いますので、基本的には各工程を上書きしています前のステップを保存します。

代わりに、ステップ数に等しいサイズのベクトルを使用するので、使用している要素を簡単に制御できます。私はデバッグの多くをやったが、私は私が間違って何をしたか理解していない最初のルーチンに使用されるもの、および第二のルーチン

ための第二1：

私は2つの異なる（しかし、基本的に同一の）機能を書きました。ここで

ソルバの2つの異なるバージョンです：

# include <iostream> 
# include <functional> 
# include <vector> 
# include <string> 
# include <cmath> 
# include <fstream> 
# include <valarray> 


using namespace std; 


std::vector<std::function<const double(const double t , std::valarray<double> u)>> func ; 


auto f2  = [](double t , double u) {return -20*u+20*sin(t)+cos(t) ; } ; 

auto numFun = [](double t , std::valarray<double> u) {return -20*u[0]+20*sin(t)+cos(t) ; } ; 


int main(){ 

     double t0 = 0.0 ; 
     double tf = 2.5 ; 
     double dt = 0.01; 


     const int N = (tf-t0)/dt; 

     func.push_back(numFun); 

     auto& f = func ; 

     double t1 ; 
     std::valarray<double> y1 ; 
     y1.resize(func.size()) ; // size = 1 
     std::valarray<double> yp1 ; 
     yp1.resize(func.size()) ; // size = 1 


     std::vector<double> t2(N+1) ; 
     std::vector<double> y2(N+1) ; 
     std::vector<double> yp2(N+1) ; 



     std::vector<double> u(1); 
     u.at(0) = 1.0 ; 

     t1 = t0 ; 

     for(auto i=0; i< u.size() ; i++) 
     { 
     y1[i]= u.at(i); 
     yp1[i] = u.at(i); 
     } 


     ofstream fn("output1.dat", ios::out); 



     for(auto i=1 ; i <=N ; i++) 
     { 

     fn << t1 << " " ; 
     for(auto j =0 ; j < y1.size() ; j++) 
      fn << y1[j] << " " ; 
     fn << endl ; 

     for(auto j=0 ; j < y1.size() ; j++) 
     { 
      yp1[j] = yp1[j] + dt * f[j](t1 , y1) ;  

      y1[j] += dt/2 * (f[j](t1 , y1) + f[j](t1+dt , yp1)); 

     } 
     t1+=dt ;  
     } 
     fn.close(); 

/* --------------------   Vector version  --------------------------------- */  

     ofstream fn2("output2.dat", ios::out); 

     t2.at(0) = t0 ; 

     y2.at(0) = u.at(0) ; 
     yp2.at(0) = u.at(0) ; 


     fn2 << t2.at(0) << " " ; 
     fn2 << y2.at(0) << " " ; 
     fn2 << endl ; 

     for(auto i=1; i <= N ; i++) 
     { 
     t2.at(i) = t2.at(i-1) + dt ; 

     // Predictor step (Exp Euler) 

     yp2.at(i) = y2.at(i-1) + dt * f2(t2.at(i-1), y2.at(i-1)) ; 

     y2.at(i) = y2.at(i-1) + dt/2 * (f2(t2.at(i-1), y2.at(i-1)) + f2(t2.at(i), yp2.at(i))) ; 


     fn2 << t2.at(i) << ' ' << y2.at(i) << " " ; 

     fn2 << endl;  
     } 

     fn2.close(); 



    return 0; 
} 

私は、出力ファイルにこの二つの異なる結果を得た：

Answer 1

このように3ベクターを使用してみてください！

  for(t=t0() ; t < tf() ; t+= dt() ) 
     { 
      f << t << ' ' ; 
      for(auto i=0 ; i< u.size() ; i++) 
       f << u[i] << " " ; 
      f << std::endl ; 

      for(auto i=0 ; i < u.size() ; i++)   
      {  
        up[i] = uc[i] + dt() * rhs.f[i](t , u) ; 


        uc[i] += dt()/2 * (rhs.f[i](t , u) + rhs.f[i](t +dt() , up)); 

        u[i] = uc[i] ; 
      }    

     } 

だから基本的には予測変数はuc[i] + ...する必要があり、最初の方法が機能するにはyp1[j] = y[j] + ...

Answer 2

の代わりに（このケースでは、彼は= yp1[j] + ...yp1[j]を書いた）標準のオイラー製剤のよう+=ことができませんでしたベクトル値のODEの場合は、3つのステップを2つ目の方法のベクトル演算に対応する3つのループに分割する必要があります。第2の方法は、y1にcopyngする前に、y1の古い値を使用して右側が最初に完全に計算されるので、y1の暗黙の中間ベクトルを含むことに注意してください。

for(auto i=1 ; i <=N ; i++) 
    { 

    fn << t1 << " " ; 
    for(auto j =0 ; j < y1.size() ; j++) 
     fn << y1[j] << " " ; 
    fn << (exp(-20*t1) + sin(t1)) << " "; 
    fn << endl ; 

    for(auto j=0 ; j < y1.size() ; j++) 
    { 
     yp1[j] = y1[j] + dt * f[j](t1 , y1) ;  

    } 
    for(auto j=0 ; j < y1.size() ; j++) 
    { 
     yc1[j] = y1[j] + dt/2 * (f[j](t1 , y1) + f[j](t1+dt , yp1)); 
    } 
    for(auto j=0 ; j < y1.size() ; j++) 
    { 
     y1[j] = yc1[j]; 
    } 
    t1+=dt ;  
    }