Intégrateurs numériques C++ pour résoudre les systèmes d'odes

Question

J'ai récemment commencé à utiliser C++ et je viens de créer une classe qui permet l'intégration d'un système d'odes défini par l'utilisateur. Il utilise deux intégrateurs différents afin de comparer ses performances. Voici la disposition générale du code:

class integrators { 
    private: 
     double ti; // initial time 
      double *xi; // initial solution 
      double tf; // end time 
      double dt; // time step 
      int  n; // number of ode's 

     public: 
      // Function prototypes 
      double f(double, double *, double *); // function to integrate 
      double rk4(int, double, double, double, double *, double *); 
      double dp8(int, double, double, double, double *, double *); 
     }; 

     // 4th Order Runge-Kutta function 
     double integrators::rk4(int n, double ti, double tf, double dt, double *xi, double *xf) { 
      // Function statements 
     } 

     // 8th Order Dormand-Prince function 
     double integrators::dp8(int n, double ti, double tf, double dt, double *xi, double *xf) { 
      // Function statements 
     } 

     // System of first order differential equations 
     double integrators::f(double t, double *x, double *dx) { 
      // Function statements 
     } 

     int main() { 
      // Initial conditions and time related parameters 
      const int n = 4; 
      double t0, tmax, dt; 
      double x0[n], xf[n]; 

      x0[0] = 0.0;       
      x0[1] = 0.0;       
      x0[2] = 1.0;    
      x0[3] = 2.0;    

      // Calling class integrators 
      integrators example01; 
      integrators example02; 

      // First integrator 
      example02.dp8(n, t0, tmax, dt, x0, xf); 

      // Second integrator 
      example01.rk4(n, t0, tmax, dt, x0, xf); 
     } 

Le problème est que le tableau contenant les conditions initiales x0 dans les principaux, les changements après l'exécution du premier intégrateur et je ne peux pas utiliser les mêmes conditions initiales pour le deuxième intégrateur, à moins que Je définis un autre tableau avec les mêmes conditions initiales (x0_rk4 et x0_dp8). Existe-t-il un moyen plus professionnel de garder ce tableau constant afin de l'utiliser dans les deux intégrateurs?

Answer 1

Non, pas vraiment. Mais il existe une solution plus élégante:

std::array<double, n> x0_rk4 = { 0.0, 0.0, 1.0, 2.0 }; 
auto x0_dp8 = x0_rk4; // copy! 

Vous devrez utiliser x0_rk4.data() pour accéder au tableau sous-jacent. Notez qu'il serait préférable que vous utilisiez std::array et d'autres fonctionnalités C++ modernes au lieu de pointeurs bruts et autres.

Answer 2

La méthode la plus simple consiste à créer une copie locale de tableau à l'intérieur des fonctions d'intégration.

Modifiez la façon dont vous transmettez 'n' à la fonction 'const int n', de sorte que vous puissiez créer quelque chose comme double currentSolution[n]; à l'intérieur et copier les éléments du tableau initial vers le nouveau. Cette approche sauvegardera votre tableau initial intact, à moins que vous ne le modifiiez "accidentellement" quelque part. Pour éviter cette probabilité de modification accidentelle, nous devons aller plus loin et utiliser un des conteneurs stl. Je pense que tout ira bien avec std::valarray<T>.

Modifiez la façon dont vous le transmettez à const std::valarray<double>& et effectuez à nouveau la copie locale non const.