Runge Kutta ODE Solver

Question

An ordinary differential equation (ODE) is an equation of the kind $$u'(x)=f(t,u(x)).$$

My program attempts to solve such ODE's numerically through explicit Runge Kutta methods. Instead of writing a new function for each and every method, it is possible to create just one function that accepts a so called butcher tableau, which contains all the necessary information for each and every Runge Kutta method.

The class

RKTableInterFace<alg,rk_Container1D, rk_Container2d>        

is such a tableau (the template parameterc rk_Container only specifies the internal data structure of the Interface and alg is the specific algorithm (e.g. euler).

If you would also like to see the rest how these classes are implemented or compile it for yourself, check the following GitHub link (at least C++11 compiler required): https://github.com/CowFreedom/ode

The algorithm follows the Wikipedia description.

Even though the code works, I have some issues:

1. Running forward in time (t_start< t_end) works but not running backwards (t_end < t_start).

2. When the end time t_end is not a multiple of the step size h it is possible that the t_end will not be evaluated.

3. I don't know if it is a good idea to use a floating point type like double for the representation of the step size h, since this sometimes leads to rounding errors when h=0.1.

4. The algorithm is just not fast. Even though you cannot test it, as it would be too much to post the header file and RKKuttaInterface.cpp, it appears to me that it would be possible to write this algorithm in a more efficient manner. Maybe parallel?

template<class fContainer, class numberContainer, class initContainer, class T, class alg, class rk_Container1D, class rk_Container2D >
auto r::Explicit_RKTemplate(r::RKTableInterface<alg, rk_Container1D,     rk_Container2D>& butcher_tableau, const fContainer& u, const numberContainer& tspan, const initContainer& u0, T h) {
    //fContainer& u contains the differential equation functions
    //numberContainer tspan contains the endpoints of the interval t_start, t_end
    //initContainer& u0 contains the initial values
    //t contains the step size

    typename numberContainer::const_iterator it_tspan = tspan.begin();
    typename fContainer::const_iterator it_fContainer = u.begin();
    typename initContainer::const_iterator it_initContainer = u0.begin();
    typename rk_Container2D::const_iterator it_rk_Table_a = butcher_tableau.getA().begin();
    typename rk_Container1D::const_iterator it_B_P = butcher_tableau.getB_P().begin();

    const initContainer::value_type t_start = *it_tspan;
    initContainer::value_type t = t_start + h;
    const initContainer::value_type t_end = *(++it_tspan);

    if (t_start > t_end) {
        h = h*-1;
    }

    std::vector<std::vector<initContainer::value_type>> trajectory_u(1, std::vector<initContainer::value_type>(u0.begin(), u0.end()));
    std::vector<initContainer::value_type> u1(trajectory_u[0].begin(), trajectory_u[0].end());
    std::vector<initContainer::value_type> trajectory_t(1, t_start);
    const size_t L = (*it_rk_Table_a).size();

    std::vector<std::vector<initContainer::value_type>> y_K_L(L, std::vector<initContainer::value_type>(u0.size(), 0));

    bool end_not_reached = true;
    double precision = std::numeric_limits<double>::denorm_min();
    std::vector<initContainer::value_type> temp(u0.size(), 0);

    rk_Container1D tk_L;

    for (size_t j = 0; j < L; j++) {
        tk_L[j] = t + (*(butcher_tableau.getC().begin() + j))*h;
    }

    for (size_t i = 1; end_not_reached; i++) {
        if (abs(t_end - t) < 2 * precision) {
            end_not_reached = false;
        }
        for (size_t s_1 = 0; s_1 < L; s_1++) {
            for (size_t s_2 = 0; s_2 <= s_1; s_2++) {
                rk_Container2D::const_iterator currentAs_1_it = butcher_tableau.getA().begin();
                for (size_t func_i = 0; func_i < u.size(); func_i++) {
                    temp[func_i] = temp[func_i] + (*((*(currentAs_1_it + s_2)).begin()))*(*std::next(it_fContainer, func_i))(y_K_L[s_1], tk_L[s_2]);
                }
            }
            for (size_t func_i = 0; func_i < u.size(); func_i++) {
                y_K_L[s_1][func_i] = trajectory_u[i - 1][func_i] + h*temp[func_i];
            }
            for (size_t func_i2 = 0; func_i2 < u.size(); func_i2++) {
                u1[func_i2] = u1[func_i2] + h*(*std::next(butcher_tableau.getB_P().begin(), s_1))*(*std::next(it_fContainer, func_i2))(y_K_L[s_1], tk_L[s_1]);
            }
            std::fill(temp.begin(), temp.end(), 0);
        }

        trajectory_u.push_back(std::vector<initContainer::value_type>(u0.size()));
        std::copy(u1.begin(), u1.end(), (trajectory_u[i]).begin());
        trajectory_t.push_back(t);
        t = t + h;
        if (t > t_end && t - h < t_end) {
            t = t - (t - t_end);
        }
    }

    return std::make_pair(trajectory_t, trajectory_u);
}

Answer 1

1. Regarding negative steps I would stop making my live hard and simply swap t_end and t_start. That way you can use the same code path for both cases without worrying about signs etc.

2. If you t_end-t_start is not a multiple of your dt you need to evaluate again with a smaller time step at the end of your computation.

3. You are stuck with float or double for h. The steps are non integral so ...

4. RK is a multistep methods, where you calculate the n-th step using the result of the n-1 th step. So it is not possible to parallelize the computation of a single step. What you can parallelize is the computation of the different ODEs.

5. You need to really step back and consider whether your code is doing the right thing and then you should ask your self how long it took you to come to that conclusion. RK is a really simple method and your code makes it just so incredible complicated

   temp[func_i] = temp[func_i] + (*((*(currentAs_1_it + s_2)).begin()))*(*std::next(it_fContainer, func_i))(y_K_L[s_1], tk_L[s_2]);
   

   There are 3 dereferences and one multiplication here and the braces are just hard to read. What makes it really bad is that it is completely unnecessary.

   *(currentAs_1_it + s_2)
   

   This is just the iterator to the s_2 elementh of the array. Why on earth do you mix index based loops with iterator loops? Either you use

   for (size_t index = 0; index < container.size(); ++index) {
       auto element = container[index];
   }  
   

   or

   for (auto it = container.begin(); it != container.end(); ++it) {
       auto element = *it;
   }  
   

   However, if you really want to simplify your code, you should use range based loops:

   for (auto&& element : container) {
   
   }  
   

   This is most likely responsible for some of the performance problems you have as you are always calculating pointers on the fly. As an example:

   for (size_t j = 0; j < L; j++) {
       tk_L[j] = t + butcher_tableau.getC()[j]*h;
   }
   

Now you see that you should also create a getC() member that take a const size_t index and returns the appropriate element.

6. I dont see where you reserve the memory for your operations. From a performance point of view this is catastrophic as you are forced to permanently reallocated huge chunks of memory

   trajectory_u.push_back(std::vector<initContainer::value_type>(u0.size()));
   std::copy(u1.begin(), u1.end(), (trajectory_u[i]).begin());
   trajectory_t.push_back(t); 
   

   I dont even know what you are trying to do here, but there is no possible reason you want to copy all that memory around. Also you know in advance the number of steps you have, but I do not see a single reserve()

So what you should do is:

   trajectory_u.reserve(numberOfSteps);
   trajectory_t.reserve(numberOfSteps);

And later:

   trajectory_u.push_back(u1);

7. You have so many typedefs, but not one for the most used ones

   std::vector<initContainer::value_type>
   

Edit

8. This code here is overly complicated. You should use short hand notations for that:

   t = t + h;
   if (t > t_end && t - h < t_end) {
       t = t - (t - t_end);
   }
   

   Should be written as:

   t += h;
   if (t > t_end) {
       t = t_end;
   }
   

   The second conditional was true by default as you increment by h

Answer 2

I have taken all your suggestions into account and updated the code! It works fabulously now!

template<class fContainer, class numberContainer, class initContainer, class T, class alg, class rk_Container1D, class rk_Container2D >

auto r::Explicit_RKTemplate(r::RKTableInterface<alg,rk_Container1D, rk_Container2D>& butcher_tableau, const fContainer& u, const numberContainer& tspan,const initContainer& u0, T h){

//fContainer& u contains the differential equation functions

//numberContainer tspan contains the endpoints of the interval t_start, t_end

//initContainer& u0 contains the initial values

//t contains the step size



typename rk_Container2D::const_iterator it_rk_Table_a = butcher_tableau.getA().begin();



const size_t u_size=u.size();

const auto t_start = *tspan.begin();

const auto t_end = *(++tspan.begin()); 



const auto t_difference=abs((t_start-t_end)/h);



std::vector<std::vector<initContainer::value_type>> trajectory_u(1,std::vector<initContainer::value_type>(u0.begin(),u0.end()));

std::vector<initContainer::value_type> trajectory_t(1,t_start);



//Check if intervall is multiple of step size to properly prereserve vector



if (floor(t_difference)==t_difference){

    trajectory_u.reserve(t_difference);

    trajectory_t.reserve(t_difference); 



}

else{

    trajectory_u.reserve(static_cast<size_t>(t_difference)+1);      

    trajectory_t.reserve(static_cast<size_t>(t_difference)+1);

}



//Array with temporary solution of next step

std::vector<initContainer::value_type> u_t(trajectory_u[0].begin(),trajectory_u[0].end());



//Number of Runge-Kutta stages

const size_t L=(*it_rk_Table_a).size();



//Solutions of intermediate locations between step sizes

std::vector<std::vector<initContainer::value_type>> y(L, std::vector<initContainer::value_type>(u_size, 0));



//Sum of the intermediate runge kutta stages

std::vector<initContainer::value_type> sum(u0.size(),0);



auto& tk_L=butcher_tableau.getC();

auto& A=butcher_tableau.getA();

auto& B=butcher_tableau.getB_P();



bool intervall_overshoot = false;

bool one_last_repetition = true;

auto t = t_start;

size_t iteration = 0;



while (t < t_end && one_last_repetition) {



    //Entscheide, ob loop nocheinmal wiederholt werden muss

    one_last_repetition = !intervall_overshoot;



    if (t+h > t_end) {

        intervall_overshoot = true;

        h = t_end - t;

    }

    //Berechne Stützstellen

    for (size_t i = 0; i < L; i++) {

        sum = u_t;

        for (size_t k = 0; k < u_size; k++) {

            for (size_t j = 0; j < i; j++) {

            //  if(k==0)

            //      std::cout << "s: " << i << "j: " << j << "Funktion:" << k <<"prior_sum[k]:"<<sum[k]<<"updated_sum[k]: "<< y[j][k] * A[i][j] <<"y[i][k]: "<< y[j][k]<< "t: "<<t<<"\n";

                sum[k] += h*y[j][k] * A[i][j];

            }

        }

        for (size_t k = 0; k < u_size; k++) {

            y[i][k] = u[k](sum, t+h*tk_L[i]);

        }

    }

    iteration++;

    t = t + h;

    //Berechne u_{t+1}

    for (size_t i = 0; i < u_size; i++) {

        for (size_t j = 0; j < L; j++) {

            u_t[i] += h*y[j][i]*B[j];

        }

    }

    trajectory_u.push_back(u_t);

    trajectory_t.push_back(t);

}

return std::make_pair(trajectory_t,trajectory_u);       

}

The entire updated code is linked under the GitHub account in the initial post.

Cheers!