Floating in C++ - part 6 - Approximating the Solution to the Cauchy Problem

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

Ordinary differential equations and Partial differential equations are two big chapters in every numerical course, so I had to spend some time on those, particularly on ODE, since there are a lot of nice and simple algorithms to implement.

Euler method

The euler method is a very simple numerical procedure for approximating the solution of a ordinary differential equation (ODE) with a given initial value.

Such problem can be described as

{ y'(t)=f(t,y(t)) y(t0)=y0

The example on Wikipedia is very descriptive, so I will show some code below:

struct explicit_euler_status{
	double t;
	double y;
	double f;
};
class explicit_euler{
	std::function<double(double, double)> f;
	explicit_euler_status st;
	double h;
public:
	explicit explicit_euler(const std::function<double(double, double)>& f_, const double y0, const double t0, double h_) : f(f_), h(h_){
		st.t = t0;
		st.y = y0;
		st.f = f(t0,y0);
	}

	explicit_euler& operator++(){
		st.y = std::fma(h, st.f, st.y);
		st.t += h;
		st.f = f(st.t, st.y);
		return *this;
	}

	const explicit_euler_status& operator*() const {
		return st;
	}
	const explicit_euler_status* operator->() const {
		return &st;
	}
};

As in other posts, I’ve implemented the algorithm with an iterator-like interface. Using this algorithm is straightforward. In the following sample, we define a Cauchy problem where we know the solution. This way we can test if our implementation really approximates the solution.

struct ode_to_solve{
	std::function<double(double,double)> f;
	std::function<double(double)> sol;
	double t0;
	double tmax;
	double y0;
};

TEST_CASE("explicit euler"){
	ode_to_solve p = {
	    [](double t, double y){return t*t *(1-3*y);},
	    [](double t){ return (1+5*std::exp(-t*t*t))/3;},
	    0,
	    3,
	    2
	};
	calcnum::explicit_euler ee(p.f, p.y0, p.t0, 0.1);
	while(ee->t < p.tmax){
		REQUIRE(Approx(p.sol(ee->t)) == ee->y);
		++ee;
	}
}

Using the same approach we can define a lot of other algorithms, for example, the Backward Euler method, Semi-implicit Euler method and others.

At one point, I noticed that the code was completely unreadable. It’s perfectly clear what a std::function<double(double,double)> is, but which one is the time parameter? The first one or the second one? Since it’s impossible to tell, and I’m too lazy to write a proper documentation, even if I would have written a wonderful one, it is surely possible to do something better. That is why I searched for another approach.

The best solution would be to give the compiler the power to stop compiling if we were using the parameters wrongly, for example if writing f(y,t) instead of f(t,y). In C++ this is very easy since we can define our own custom types:

struct cauchy_time{
	double t;
}

If we change the function signature to std::function<double(cauchy_time,double)>, there is no need to document the difference between the parameters, and the compiler will warn us if we are interchanging them. Unfortunately doing so means that we need to do an explicit conversion for every function that needs to take our cauchy_time as a parameter but accepts a double, unless we will rewrite those functions and all operator like operator+, operator+=, operator==.

Therefore, I decided to use a simple alias: using cauchy_time = double;. This means that double and cauchy_time have exactly the same type and the compiler will not issue a compile error if we write f(y,t) instead of f(t,y). The upside is that we can write std::function<double(cauchy_time,double)> where it is relevant, and the code will be much more self-documenting than before.

If you think that such a struct or alias could be totally avoided, you are right, but consider following: for the [Backward Euler] method we need a callback function:

std::function<double(double, double, std::function<double(double,double)>, double)> callback;

I have no idea what double stands for in the different positions, and I’m the one who wrote it! I bet no one will ever be able to guess it. Now consider

using cauchy_time = double;
using cauchy_var = double;
using cauchy_eq = std::function<cauchy_var(cauchy_time,cauchy_var)>;
std::function<cauchy_var(cauchy_var, cauchy_time, cauchy_eq, double)> callback;

It may still not be obvious how such a function should work, but you have an idea of the semantical meaning of the parameters and how they can possibly interact.

As a bonus point, the compiler can generate easier to understand error messages, then the ones generate if using only the double type and awesome documentation:

candidate constructor not viable: no known conversion from 'const double' to 'std::function<cauchy_var (cauchy_var, cauchy_time, cauchy_eq, double)>' (aka 'function<double (double, double, function<double (double, double)>, double)>') for 3rd argument

The real question is, why would someone need such a definition? The downside of an implicit method is that at every iteration, the solution needs to be approximated by a second method. If we do not want to hardcode the second method, but let the caller provide it (at runtime or compile time, it does not matter), we need something like a callback. We can combine, for example, the fixed point iteration method with the implicit Euler method construct in order to solve an ODE:

cauchy_var my_callback(cauchy_var yn, cauchy_time tn1, cauchy_eq f, double h){
	auto fun = [yn, tn1, f, h](cauchy_var y){
		return h* f(tn1,y) + yn;
	};
	calcnum::fixed_point fp(fun, yn);
	for(int i = 0; i != 100; ++i){
		++fp;
	}
	return fp->f;
}

TEST_CASE("implicit euler"){
	cauchy_prob cp = {
	    [](double t, double y){return t*t*(1-3*y);},
	    cauchy_time(0),
	    cauchy_var(2)
	};
	double tmax = 3;
	auto sol = [](double t){ return (1+5*std::exp(-t*t*t))/3;};

	const double h = 0.01;
	calcnum::implicit_euler ee(cp, h, my_callback);
	double err = 0;
	while(ee->t < tmax){
		REQUIRE(sol(ee->t) == Approx(ee->y).epsilon(h));
		++ee;
	}
}

this way, at every iteration, the function my_callback is automatically called for solving the implicit part of the algorithm.