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(t0) = y0
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(t0) = T0 cools down in an environment of temperature TR:
or
It says that the cooling rate 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
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 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.
#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.