user
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 {\displaystyle {\frac {dy(t)}{dt}}=f(t,y(t))}}

with an initial value

:   :   :   y
            (

            t

            0
            )
            =

            y

            0
            {\displaystyle {\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 {\displaystyle {\frac {dy(t)}{dt}}\approx {\frac {y(t+h)-y(t)}{h}}}}

then solve for 

y
(
t
+
h
)
{\displaystyle {\displaystyle y(t+h)}}
:

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

            d
            y
            (
            t
            )

            d
            t
            {\displaystyle {\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 {\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 {\displaystyle y\_{n+1}=y\_{n}+h\,f(t\_{n},y\_{n})}}

where   

h
{\displaystyle {\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 {\displaystyle T(t\_{0})=T\_{0}}}
   cools down in an environment of temperature   

T

R
{\displaystyle {\displaystyle T\_{R}}}
:

:   :   :   d
            T
            (
            t
            )

            d
            t
            =
            −
            k

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

or

:   :   :   d
            T
            (
            t
            )

            d
            t
            =
            −
            k

            (
            T
            (
            t
            )
            −

            T

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

It says that the cooling rate   

d
T
(
t
)

d
t
{\displaystyle {\displaystyle {\frac {dT(t)}{dt}}}}
   of the object is proportional to the current temperature difference   

Δ
T
=
(
T
(
t
)
−

T

R
)
{\displaystyle {\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 {\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 {\displaystyle T\_{0}}}
                 shall be   100 °C
            * room temperature   

              T

              R
              {\displaystyle {\displaystyle T\_{R}}}
                 shall be   20 °C
            * cooling constant     

              k
              {\displaystyle {\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.
assistant
(formerly Perl 6)

```
sub euler ( &f, $y0, $a, $b, $h ) {
    my $y = $y0;
    my @t_y;
    for $a, * + $h ... * > $b -> $t {
        @t_y[$t] = $y;
        $y += $h × f( $t, $y );
    }
    @t_y
}

constant COOLING_RATE = 0.07;
constant AMBIENT_TEMP =   20;
constant INITIAL_TEMP =  100;
constant INITIAL_TIME =    0;
constant FINAL_TIME   =  100;

sub f ( $time, $temp ) {
    -COOLING_RATE × ( $temp - AMBIENT_TEMP )
}

my @e;
@e[$_] = euler( &f, INITIAL_TEMP, INITIAL_TIME, FINAL_TIME, $_ ) for 2, 5, 10;

say 'Time Analytic   Step2   Step5  Step10     Err2     Err5    Err10';

for INITIAL_TIME, * + 10 ... * >= FINAL_TIME -> $t {

    my $exact = AMBIENT_TEMP + (INITIAL_TEMP - AMBIENT_TEMP)
                              × (-COOLING_RATE × $t).exp;

    my $err = sub { @^a.map: { 100 × ($_ - $exact).abs / $exact } }

    my ( $a, $b, $c ) = map { @e[$_][$t] }, 2, 5, 10;

    say $t.fmt('%4d '), ( $exact, $a, $b, $c )».fmt(' %7.3f'),
                            $err([$a, $b, $c])».fmt(' %7.3f%%');
}

```

Output:

```
Time Analytic   Step2   Step5  Step10     Err2     Err5    Err10
   0  100.000 100.000 100.000 100.000   0.000%   0.000%   0.000%
  10   59.727  57.634  53.800  44.000   3.504%   9.923%  26.331%
  20   39.728  37.704  34.281  27.200   5.094%  13.711%  31.534%
  30   29.797  28.328  26.034  22.160   4.927%  12.629%  25.629%
  40   24.865  23.918  22.549  20.648   3.808%   9.313%  16.959%
  50   22.416  21.843  21.077  20.194   2.555%   5.972%   9.910%
  60   21.200  20.867  20.455  20.058   1.569%   3.512%   5.384%
  70   20.596  20.408  20.192  20.017   0.912%   1.959%   2.808%
  80   20.296  20.192  20.081  20.005   0.512%   1.057%   1.432%
  90   20.147  20.090  20.034  20.002   0.281%   0.559%   0.721%
 100   20.073  20.042  20.014  20.000   0.152%   0.291%   0.361%
```