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:

\frac{dy(t)}{dt} = f(t,y(t))

with an initial value

y(t0) = y0

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

\frac{dy(t)}{dt} \approx \frac{y(t+h)-y(t)}{h}

then solve for y(t + h):

y(t+h) \approx y(t) + h \, \frac{dy(t)}{dt}

which is the same as

y(t+h) \approx y(t) + h \, f(t,y(t))

The iterative solution rule is then:

y_{n+1} = y_n + h \, f(t_n, y_n)

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(t0) = T0 cools down in an environment of temperature TR:

\frac{dT(t)}{dt} = -k \, \Delta T


\frac{dT(t)}{dt} = -k \, (T(t) - T_R)

It says that the cooling rate \frac{dT(t)}{dt} of the object is proportional to the current temperature difference ΔT = (T(t) − TR) to the surrounding environment.

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

T(t) = T_R + (T_0 - T_R) \; e^{-k t}


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 T0 shall be 100 °C, the room temperature TR 20 °C, and the cooling constant k 0.07. The time interval to calculate shall be from 0 s to 100 s.

Euler Method Newton Cooling.png

#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


Content is available under GNU Free Documentation License 1.2.