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.


program euler_method
use iso_fortran_env, only: real64
implicit none
 
abstract interface
  ! a derivative dy/dt as function of y and t
  function derivative(y, t)
    use iso_fortran_env, only: real64
    real(real64) :: derivative
    real(real64), intent(in) :: t, y
  end function
end interface
 
real(real64), parameter :: T_0 = 100, T_room = 20, k = 0.07, a = 0, b = 100, &
    h(3) = [2.0, 5.0, 10.0]
 
integer :: i
 
! loop over all step sizes
do i = 1, 3
  call euler(newton_cooling, T_0, a, b, h(i))
end do
 
contains
 
! Approximates y(t) in y'(t) = f(y, t) with y(a) = y0 and t = a..b and the 
! step size h.
subroutine euler(f, y0, a, b, h)
  procedure(derivative) :: f
  real(real64), intent(in) :: y0, a, b, h
  real(real64) :: t, y
 
  if (a > b) return
  if (h <= 0) stop "negative step size"
 
  print '("# h = ", F0.3)', h
 
  y = y0
  t = a
 
  do
    print *, t, y
    t = t + h
    if (t > b) return
    y = y + h * f(y, t)
  end do
end subroutine
 
 
! Example: Newton's cooling law, f(T, _) = -k*(T - T_room) 
function newton_cooling(T, unused) result(dTdt)
  real(real64) :: dTdt
  real(real64), intent(in) :: T, unused
  dTdt = -k * (T - T_room)
end function
 
end program