C++: Euler Method

Euler’s method numerically approximates solutions of first-order ordinary differential equations (ODEs) with a given initial value. It is an explicit method for solving initial value problems (IVPs), as described in the wikipedia page. The ODE has to be provided in the following form:

with an initial value

y(t₀) = y₀

To get a numeric solution, we replace the derivative on the LHS with a finite difference approximation:

then solve for y(t + h):

which is the same as

The iterative solution rule is then:

h is the step size, the most relevant parameter for accuracy of the solution. A smaller step size increases accuracy but also the computation cost, so it has always has to be hand-picked according to the problem at hand.

Example: Newton’s Cooling Law

Newton’s cooling law describes how an object of initial temperature T(t₀) = T₀ cools down in an environment of temperature T_R:

or

It says that the cooling rate of the object is proportional to the current temperature difference ΔT = (T(t) − T_R) to the surrounding environment.

The analytical solution, which we will compare to the numerical approximation, is

Task

The task is to implement a routine of Euler’s method and then to use it to solve the given example of Newton’s cooling law with it for three different step sizes of 2 s, 5 s and 10 s and to compare with the analytical solution. The initial temperature T₀ shall be 100 °C, the room temperature T_R 20 °C, and the cooling constant k 0.07. The time interval to calculate shall be from 0 s to 100 s.

#include <iomanip>
#include <iostream>

typedef double F(double,double);

/*
Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and
t=a..b and the step size h.
*/
void euler(F f, double y0, double a, double b, double h)
{
	double y = y0;
	for (double t = a; t < b; t += h)
	{
		std::cout << std::fixed << std::setprecision(3) << t << " " << y << "\n";
		y += h * f(t, y);
	}
	std::cout << "done\n";
}

// Example: Newton's cooling law
double newtonCoolingLaw(double, double t)
{
	return -0.07 * (t - 20);
}

int main()
{
	euler(newtonCoolingLaw, 100, 0, 100,  2);
	euler(newtonCoolingLaw, 100, 0, 100,  5);
	euler(newtonCoolingLaw, 100, 0, 100, 10);
}

Last part of output:

...
0.000 100.000
10.000 44.000
20.000 27.200
30.000 22.160
40.000 20.648
50.000 20.194
60.000 20.058
70.000 20.017
80.000 20.005
90.000 20.002
done

SOURCE

Content is available under GNU Free Documentation License 1.2.