C++ numerical integrators to solve systems of ode&

Question

I recently started using C++ and I just created a class that allows the integration of a user-defined system of ode's. It uses two different integrators in order to compare its performance. Here is the general layout of the 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);
        }

The problem is that the array containing the initial conditions x0 in main, changes after executing the first integrator and I cannot use the same initial conditions for the second integrator, unless I define another array with the same initial conditions (x0_rk4 and x0_dp8). Is there a more professional way to keep this array constant in order to use it in both integrators?

Answer 1

No, not really. But there exist a more elegant solution:

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

You will have to use x0_rk4.data() to access the underlying array. Note that it would be better if you used std::array and other modern C++ features instead of raw pointers and the like.

Answer 2

The easiest way is to make a local copy of array inside integrating functions.

Change the way you are passing 'n' to function to 'const int n', so you can make something like double currentSolution[n]; inside and copy elements from initial array to new one. This approach will save your initial array intact, unless you will "accidently" modify it somewhere.

To prevent this probability of accident modification we need to go deeper and use one of stl containers. I think you will be fine with std::valarray<T>.

Change the way you are passing it to const std::valarray<double>& and again make non-const local copy.