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:

 






d
y
(
t
)


d
t



=
f
(
t
,
y
(
t
)
)


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


with an initial value

 



y
(

t

0


)
=

y

0




{\displaystyle y(t_{0})=y_{0}}


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

 






d
y
(
t
)


d
t



≈



y
(
t
+
h
)
−
y
(
t
)

h




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


then solve for 



y
(
t
+
h
)


{\displaystyle y(t+h)}

:

 



y
(
t
+
h
)
≈
y
(
t
)
+
h




d
y
(
t
)


d
t





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


which is the same as

 



y
(
t
+
h
)
≈
y
(
t
)
+
h

f
(
t
,
y
(
t
)
)


{\displaystyle 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


)


{\displaystyle y_{n+1}=y_{n}+h\,f(t_{n},y_{n})}


where   



h


{\displaystyle 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




{\displaystyle T(t_{0})=T_{0}}

   cools down in an environment of temperature   




T

R




{\displaystyle T_{R}}

:

 






d
T
(
t
)


d
t



=
−
k

Δ
T


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


or

 






d
T
(
t
)


d
t



=
−
k

(
T
(
t
)
−

T

R


)


{\displaystyle {\frac {dT(t)}{dt}}=-k\,(T(t)-T_{R})}




It says that the cooling rate   






d
T
(
t
)


d
t





{\displaystyle {\frac {dT(t)}{dt}}}

   of the object is proportional to the current temperature difference   



Δ
T
=
(
T
(
t
)
−

T

R


)


{\displaystyle \Delta T=(T(t)-T_{R})}

   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




{\displaystyle T(t)=T_{R}+(T_{0}-T_{R})\;e^{-kt}}




Task

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.



Initial values

   initial temperature   




T

0




{\displaystyle T_{0}}

   shall be   100 °C
   room temperature   




T

R




{\displaystyle T_{R}}

   shall be   20 °C
   cooling constant     



k


{\displaystyle k}

     shall be   0.07 
   time interval to calculate shall be from   0 s   ──►   100 s


A reference solution (Common Lisp) can be seen below.   We see that bigger step sizes lead to reduced approximation accuracy.


#include <stdio.h>
#include <math.h>
 
typedef double (*deriv_f)(double, double);
#define FMT " %7.3f"
 
void ivp_euler(deriv_f f, double y, int step, int end_t)
{
	int t = 0;
 
	printf(" Step %2d: ", (int)step);
	do {
		if (t % 10 == 0) printf(FMT, y);
		y += step * f(t, y);
	} while ((t += step) <= end_t);
	printf("\n");
}
 
void analytic()
{
	double t;
	printf("    Time: ");
	for (t = 0; t <= 100; t += 10) printf(" %7g", t);
	printf("\nAnalytic: ");
 
	for (t = 0; t <= 100; t += 10)
		printf(FMT, 20 + 80 * exp(-0.07 * t));
	printf("\n");
}
 
double cooling(double t, double temp)
{
	return -0.07 * (temp - 20);
}
 
int main()
{
	analytic();
	ivp_euler(cooling, 100, 2, 100);
	ivp_euler(cooling, 100, 5, 100);
	ivp_euler(cooling, 100, 10, 100);
 
	return 0;
}