Search code examples
c++differential-equations

what is the difference in this 2 code?

I wrote two codes for solving a differential problem, both using predictor corrector approach.

The first one is for solving differential system as well, so basically the it uses a vector y of size equal to the number of the equation to solve in the system, so this one basically is overwritten each step and doesn't store the step before.

The second one instead uses a vector of size equal to the number of steps, so it's easy to control which elements we are using.

I wrote two different (but basically identical) functions: one used for the first routine, and the second one for the second routine

I did a lot of debug but I don't understand what I did wrong.

Here are the two different versions of the solver:

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

I got this two different result in the output file : this

Solution

  • try to use 3 vector in this way ! 

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

         }

    So basically the predictor variable have to be uc[i] + ... and could not be += like in the standard Euler formulation (in this case he wrote yp1[j] = yp1[j] + ... ) instead of yp1[j] = y[j] + ...