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*_{0}) =*y*_{0}

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*_{0}) = *T*_{0} 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*_{0} 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

Content is available under GNU Free Documentation License 1.2.