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:
Running forward in time (t_start< t_end) works but not running backwards (t_end < t_start).
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.
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.
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);
}