Floating in C++ - part 7 - Runge-Kutta

The introduction to this article series can be found here, whereas here you can find all posts.

In the last article I mentioned how there are many algorithms for solving ODE’s, but showed only the implementation for explicit and implicit Euler. The reason behind my decision was that those algorithms, and many others as well, are part of a larger family of algorithm: The Runge-Kutta methods.

Again, Wikipedia does a great job by describing how to put up such a system. TLDR: The Runga-Kutta methods approximate the cauchy problem through a set of equations. By putting those equations together we obtain the Runge-Kutta Matrix, and thank to that, approximating the Cauchy problem is equivalent to solve a linear system.

The equations of the Runge-Kutta methods are uniquely identified by a set of coefficients, weights, and nodes. Those are commonly grouped together in a so-called Butcher tableau.

I’ve collected some tableaux found in textbooks and online. A nice list can also be found here. The most famous tableaux are surely Forward Euler, Backward Euler and the midpoint method, since they are very simple to explain and implement by hand.

Explicit Methods

The Runga-Kutta methods can be splitted in three main families: explicit, semi-implicit and implicit. Then there are adaptive and non-adaptive methods…​ in fact one can put together different approximations techniques, prove that they converge to the solutions, and rediscover a new family of approximation methods.

The easiest family of method to implement is the explicit method. It consists of a triangular matrix system, therefore no need for a particularly robust matrix system for solving those equations. I did not implement the semi-implicit and implicit methods as I believe that a third party library with a robust matrix system is required for getting good results.

If we are only interested in one-dimensional systems, the code will look like

std::array<double, N> K = {};
for(std::size_t i = 0; i != N; ++i){
	double kb = 0;
	for(std::size_t j = 0; j != i; ++j){ // avoid j>i since K[j] == 0
		kb += K[i] * B[{i,j}];
	}
	double t_i = std::fma(a[i], h, t);
	K[i] = f(t_i, y + h*kb);
}

double kc = 0;
for(std::size_t i = 0; i != N; ++i){
	kc += K[i] * c[i];
}
y += h*kc;
t += h;

whereas for a multidimensional system, the code will look similar to this

calcnum::matrix<yN, N> K;
for(std::size_t i = 0; i != N; ++i){
	double t_i = std::fma(a[i], h, t);
	auto Ki = f(t_i, y + h*K*Bt[i] );
	set_col(K, i, Ki);
}
y += h*K*c;
t += h;

The code looks extremely compact and similar to the equations that can be found on most textbooks. That is thanks to the language that allows to provide custom multiplication and sum operators, auto, function objects, function and operator overloading. If we strip those features from C++ (aka go back to C), we would need to write something like

calcnum::matrix<yN, N> K;
for(std::size_t i = 0; i != N; ++i){
	double t_i = std::fma(get_elemen(a, i), h, t);
	calcnum::vector<yN> hKB_i = mult_mat_by_scalar(h, mult_mat_by_vec(K, get_row_vec(Bt, i), yN, N), yN);
	calcnum::vector<yN> y_i = sum_vec(y, hKB_im yN);
	calcnum::vector<yN> Ki = apply_fun(f, t_i, y_i, yN);
	set_col(K, i, Ki, yN);
}
y = sum_vec(y, mult_mat_by_scalar(h, mult_mat_by_vec(K,c,yN), yN, N), yN);
t += h;

In fact, in C it’s a common convention to pass pointers to matrices and vectors, and not return them by value. So the code would resemble even less the mathematical formulae. And its always confusing to remember which parameter stands for the number of rows, and which stands for the number of columns. The type system of C++ makes it so much easier to avoid these kind errors, but unfortunately dynamic languages like Octave, that provide a very good syntax, can not take advantage of this information. The information about the sizes of matrices is only used and checked when the interpreter is running the required expression.

There seems to be a linter for Matlab, which I’ve never tried it. From what can be gathered from the examples it does not seem able to verify the correctness of the code, but it provides suggestions on how to improve it.

I guess that a linter that checks for matrix and vector sizes could be possible, even if it might not get every use case, it could avoid many common errors.

Going back to {cpp] and numeric calculus: My implementation can be found on GitHub: I left both the one-dimensional as the multidimensional explicit Runge-Kutta iterator.

Adaptive Methods

The Runga-Kutta methods can also be categorized by other criteria, and not only between explicit, semi-implicit and implicit.

An important set of methods are the adaptive ones. Adaptive methods are iterative methods that analyze how good hey are approximating the problem, and decide to take smaller or larger steps.

This is important in the case of stiff equations. Unless the step size is very small, the approximated solution will be completely different from the exact solution. Since we would like our algorithms to finish, using a minuscule step is not a viable solution, but also using a step size that is too big would not be helpful either.

A first naive idea is to calculate the solution with two different step sizes, for example h and 0.5*h. On one hand, if the solutions are similar, then the equation is probably not stiff, and we can use a larger step. On the other hand, if the solutions are not similar, we should probably use a smaller step.

The implementation can be found, like all other algorithms, on GitHub.

Fehlberg Adaptive Methods

One of the most expensive part of the Runge-Kutta method is calculating the K matrix, and, without further assumptions, for every step this matrix is different. If we give a closer look the the equations, or at our implementations, we may notice that it does not depend on the weights at all. This means that we could reuse the same Matrix if we could identify two Butcher tableaux which differ only in the weights. And that’s the idea behind the Runge–Kutta–Fehlberg method.

Compared to the previous adaptive methods, the algorithm will have a better performance because we will need roughly a third of matrix computations.

What is a good solution for solving the same problem with different weights? Supposing that both methods give us an approximation of the solution, we have information about the actual error. By using that said information, it is possible to have a better control on the complete error, by resizing the size of the steps.

Since the two methods have a different estimation order, when we have a difference under a certain threshold, we can use a larger step. If not, we should probably use a smaller step.

while(t < tmax){
	auto err = std::numeric_limits<double>::infinity();
	for(int i = 0; i != nitmax; ++i){
		const auto hmin2 = std::min(hmin, tmax-t);
		h = (i == nitmax-1) ? hmin2 : calcnum::fclamp(h*q, hmin2, hmax);

		cauchy_probN<1> cpN_sub = {
		    cpN.f,
		    t,
		    y
		};

		calcnum::rk_exp_iter<rk_size> rk(cpN_sub, fehl_b.a, Bt, fehl_b.c1, h);
		++rk;

		const auto hKc_diff = h *rk->K*cdiff;
		err = calcnum::norm_inf(hKc_diff);
		q = std::pow( toll/err , 1/d(fehl_b.order+1) );

		if(err <= toll){ // set y, t and exit inner loop
			y = rk->y + hKc_diff;
			t = rk->t;
			break;
		}
	}
}

As shown in the example, at every step a solution to the Runge-Kutta system gets approximated. Through the matrix K and the weights it is possible to have an approximation of the local error between two methods.

Compared to the first adaptive method

while(t < tmax){
	auto err = std::numeric_limits<double>::infinity();
	for(int i = 0; i != nitmax; ++i){
		const auto hmin2 = std::min(hmin, tmax-t); // adapt hmin to avoid overflow
		REQUIRE(hmin2 > 0);
		h = (i == nitmax-1) ? hmin2 : calcnum::fclamp(h*q, hmin2, hmax);
		cauchy_probN<1> cpN_sub = {
		    cpN.f,
		    t,
		    y
		};

		calcnum::rk_exp_iter<rk_size> rk(cpN_sub, fehl_b.a, Bt, fehl_b.c1, h);
		calcnum::rk_exp_iter<rk_size> rk2(cpN_sub, fehl_b.a, Bt, fehl_b.c1, h/2);
		++rk;
		++rk2; ++rk2;

		err = calcnum::norm_inf(rk->y - rk2->y)/(std::pow(2.0,p)-1);
		q = std::pow(toll/(err*std::pow(2,p)), 1/d(p+1));

		if(err <= toll){ // set y, t and exit inner loop
			y = rk2->y;
			t = rk->t;
			break;
		}
	}
}

We are doing only one third of the work. It might not look like a big difference, but at every rk` and `rk2, we are solving a matrix system of size rk_size x rk_size. And we are solving those systems at every step at least(!) one time in order to see if we need to change our step size.