Euler method

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:

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

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:

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

then solve for $y(t+h)$ :

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

which is the same as

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})

$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}$ :

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

or

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

It says that the cooling rate ${\frac {dT(t)}{dt}}$ of the object is proportional to the current temperature difference $\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^{-kt}

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.

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

Ada

The solution is generic, usable for any floating point type. The package specification: <lang Ada> generic

  type Number is digits <>;

package Euler is

  type Waveform is array (Integer range <>) of Number;
  function Solve
           (  F      : not null access function (T, Y : Number) return Number;
              Y0     : Number;
              T0, T1 : Number;
              N      : Positive
           )  return Waveform;

end Euler; </lang> The function Solve returns the solution of the differential equation for each of N+1 points, starting from the point T0. The implementation: <lang Ada> package body Euler is

  function Solve
           (  F      : not null access function (T, Y : Number) return Number;
              Y0     : Number;
              T0, T1 : Number;
              N      : Positive
           )  return Waveform is
     dT : constant Number := (T1 - T0) / Number (N);
  begin
     return Y : Waveform (0..N) do
        Y (0) := Y0;
        for I in 1..Y'Last loop
           Y (I) := Y (I - 1) + dT * F (T0 + dT * Number (I - 1), Y (I - 1));
        end loop;
     end return;
  end Solve;

end Euler; </lang> The test program: <lang Ada> with Ada.Text_IO; use Ada.Text_IO; with Euler;

procedure Test_Euler_Method is

  package Float_Euler is new Euler (Float);
  use Float_Euler;

  function Newton_Cooling_Law (T, Y : Float) return Float is
  begin
     return -0.07 * (Y - 20.0);
  end Newton_Cooling_Law;
  
  Y : Waveform := Solve (Newton_Cooling_Law'Access, 100.0, 0.0, 100.0, 10);

begin

  for I in Y'Range loop
     Put_Line (Integer'Image (10 * I) & ":" & Float'Image (Y (I)));
  end loop;

end Test_Euler_Method; </lang> Sample output:

 0: 1.00000E+02
 10: 4.40000E+01
 20: 2.72000E+01
 30: 2.21600E+01
 40: 2.06480E+01
 50: 2.01944E+01
 60: 2.00583E+01
 70: 2.00175E+01
 80: 2.00052E+01
 90: 2.00016E+01
 100: 2.00005E+01

ALGOL 68

Translation of: D

Note: This specimen retains the original D coding style.

Works with: ALGOL 68 version Revision 1 - no extensions to language used.

Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny.

<lang algol68># Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and t=a..b and the step size h.

PROC euler = (PROC(REAL,REAL)REAL f, REAL y0, a, b, h)REAL: (

   REAL y := y0,
        t := a;
   WHILE t < b DO
     printf(($g(-6,3)": "g(-7,3)l$, t, y));
     y +:= h * f(t, y);
     t +:= h
   OD;
   printf($"done"l$);
   y

);

Example: Newton's cooling law #

PROC newton cooling law = (REAL time, t)REAL: (

   -0.07 * (t - 20)

);

main: (

  euler(newton cooling law, 100, 0, 100,  10)

)</lang> Ouput:

 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

C

<lang c>#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; }</lang>output<lang> Time: 0 10 20 30 40 50 60 70 80 90 100 Analytic: 100.000 59.727 39.728 29.797 24.865 22.416 21.200 20.596 20.296 20.147 20.073

Step  2:  100.000  57.634  37.704  28.328  23.918  21.843  20.867  20.408  20.192  20.090  20.042
Step  5:  100.000  53.800  34.280  26.034  22.549  21.077  20.455  20.192  20.081  20.034  20.014
Step 10:  100.000  44.000  27.200  22.160  20.648  20.194  20.058  20.017  20.005  20.002  20.000</lang>

C++

Translation of: D

<lang cpp>#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);

}</lang> 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

C#

<lang csharp>using System;

namespace prog { class MainClass { const float T0 = 100f; const float TR = 20f; const float k = 0.07f; readonly static float[] delta_t = {2.0f,5.0f,10.0f}; const int n = 100;

public delegate float func(float t); static float NewtonCooling(float t) { return -k * (t-TR); }

public static void Main (string[] args) { func f = new func(NewtonCooling); for(int i=0; i<delta_t.Length; i++) { Console.WriteLine("delta_t = " + delta_t[i]); Euler(f,T0,n,delta_t[i]); } }

public static void Euler(func f, float y, int n, float h) { for(float x=0; x<=n; x+=h) { Console.WriteLine("\t" + x + "\t" + y); y += h * f(y); } } } }</lang>

Clay

<lang Clay> import printer.formatter as pf;

euler(f, y, a, b, h) {

   while (a < b) {
       println(pf.rightAligned(2, a), " ", y);
       a += h;
       y += h * f(y);
   }

}

main() {

   for (i in [2.0, 5.0, 10.0]) {
       println("\nFor delta = ", i, ":");
       euler((temp) => -0.07 * (temp - 20), 100.0, 0.0, 100.0, i);
   }

} </lang>

Example output:

For delta = 10:
 0 100
10 43.99999999999999
20 27.2
30 22.16
40 20.648
50 20.1944
60 20.05832
70 20.017496
80 20.0052488
90 20.00157464

Common Lisp

<lang lisp>;; 't' usually means "true" in CL, but we need 't' here for time/temperature. (defconstant true 'cl:t) (shadow 't)

Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and t=a..b and the step size h.

(defun euler (f y0 a b h)

 ;; Set the initial values and increments of the iteration variables.
 (do ((t a  (incf t h))
      (y y0 (incf y (* h (funcall f t y)))))

     ;; End the iteration when t reaches the end b of the time interval.
     ((>= t b) 'DONE)

     ;; Print t and y(t) at every step of the do loop.
     (format true "~6,3F  ~6,3F~%" t y)))

Example

Newton's cooling law, f(t,T) = -0.07*(T-20)

(defun newton-cooling (time T) (* -0.07 (- T 20)))

Generate the data for all three step sizes (2,5 and 10).

(euler #'newton-cooling 100 0 100 2) (euler #'newton-cooling 100 0 100 5) (euler #'newton-cooling 100 0 100 10)</lang>

Example 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

D

<lang d>import std.stdio, std.range;

/** 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)(in F f, in double y0,

             in double a, in double b, in double h) {
   double y = y0;
   foreach (t; iota(a, b, h)) {
       writefln("%.3f  %.3f", t, y);
       y += h * f(t, y);
   }
   writeln("done");

}

void main() {

   /// Example: Newton's cooling law
   static newtonCoolingLaw(in double time, in double t) {
       return -0.07 * (t - 20);
   }

   euler(&newtonCoolingLaw, 100, 0, 100,  2);
   euler(&newtonCoolingLaw, 100, 0, 100,  5);
   euler(&newtonCoolingLaw, 100, 0, 100, 10);

}</lang> Last part of the 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

Fortran

Works with: Fortran version 2008

<lang fortran>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</lang> Output for h = 10:

# h = 10.000
   0.0000000000000000        100.00000000000000     
   10.000000000000000        43.999999761581421     
   20.000000000000000        27.199999856948853     
   30.000000000000000        22.159999935626985     
   40.000000000000000        20.647999974250794     
   50.000000000000000        20.194399990344049     
   60.000000000000000        20.058319996523856     
   70.000000000000000        20.017495998783350     
   80.000000000000000        20.005248799582862     
   90.000000000000000        20.001574639859214     
   100.00000000000000        20.000472391953071

Go

<lang go>package main

import (

   "fmt"
   "math"

)

// fdy is a type for function f used in Euler's method. type fdy func(float64, float64) float64

// eulerStep computes a single new value using Euler's method. // Note that step size h is a parameter, so a variable step size // could be used. func eulerStep(f fdy, x, y, h float64) float64 {

   return y + h*f(x, y)

}

// Definition of cooling rate. Note that this has general utility and // is not specific to use in Euler's method.

// newCoolingRate returns a function that computes cooling rate // for a given cooling rate constant k. func newCoolingRate(k float64) func(float64) float64 {

   return func(deltaTemp float64) float64 {
       return -k * deltaTemp
   }

}

// newTempFunc returns a function that computes the analytical solution // of cooling rate integrated over time. func newTempFunc(k, ambientTemp, initialTemp float64) func(float64) float64 {

   return func(time float64) float64 {
       return ambientTemp + (initialTemp-ambientTemp)*math.Exp(-k*time)
   }

}

// newCoolingRateDy returns a function of the kind needed for Euler's method. // That is, a function representing dy(x, y(x)). // // Parameters to newCoolingRateDy are cooling constant k and ambient // temperature. func newCoolingRateDy(k, ambientTemp float64) fdy {

   crf := newCoolingRate(k)
   // note that result is dependent only on the object temperature.
   // there are no additional dependencies on time, so the x parameter
   // provided by eulerStep is unused.
   return func(_, objectTemp float64) float64 {
       return crf(objectTemp - ambientTemp)
   }

}

func main() {

   k := .07
   tempRoom := 20.
   tempObject := 100.
   fcr := newCoolingRateDy(k, tempRoom)
   analytic := newTempFunc(k, tempRoom, tempObject)
   for _, deltaTime := range []float64{2, 5, 10} {
       fmt.Printf("Step size = %.1f\n", deltaTime)
       fmt.Println(" Time Euler's Analytic")
       temp := tempObject
       for time := 0.; time <= 100; time += deltaTime {
           fmt.Printf("%5.1f %7.3f %7.3f\n", time, temp, analytic(time))
           temp = eulerStep(fcr, time, temp, deltaTime)
       }
       fmt.Println()
   }

}</lang> Output, truncated:

...
 85.0  20.053  20.208
 90.0  20.034  20.147
 95.0  20.022  20.104
100.0  20.014  20.073

Step size = 10.0
 Time Euler's Analytic
  0.0 100.000 100.000
 10.0  44.000  59.727
 20.0  27.200  39.728
 30.0  22.160  29.797
 40.0  20.648  24.865
 50.0  20.194  22.416
 60.0  20.058  21.200
 70.0  20.017  20.596
 80.0  20.005  20.296
 90.0  20.002  20.147
100.0  20.000  20.073

Groovy

Generic Euler Method Solution

The following is a general solution for using the Euler method to produce a finite discrete sequence of points approximating the ODE solution for y as a function of x.

In the eulerStep closure argument list: x_n and y_n together are the previous point in the sequence. h is the step distance to the next point's x value. dydx is a closure representing the derivative of y as a function of x, itself expressed (as required by the method) as a function of x and y(x).

The eulerMapping closure produces an entire approximating sequence, expressed as a Map object. Here, x₀ and y₀ together are the first point in the sequence, the ODE initial conditions. h and dydx are again the step distance and the derivative closure. stopCond is a closure representing a "stop condition" that causes the the eulerMapping closure to stop after a finite number of steps; the given default value causes eulerMapping to stop after 100 steps. <lang groovy>def eulerStep = { xn, yn, h, dydx ->

   (yn + h * dydx(xn, yn)) as BigDecimal

}

Map eulerMapping = { x0, y0, h, dydx, stopCond = { xx, yy, hh, xx0 -> abs(xx - xx0) > (hh * 100) }.rcurry(h, x0) ->

   Map yMap = [:]
   yMap[x0] = y0 as BigDecimal
   def x = x0
   while (!stopCond(x, yMap[x])) {
       yMap[x + h] = eulerStep(x, yMap[x], h, dydx)
       x += h
   }
   yMap

} assert eulerMapping.maximumNumberOfParameters == 5</lang>

Specific Euler Method Solution for the "Temperature Diffusion" Problem (with Newton's derivative formula and constants for environment temperature and object conductivity given) <lang groovy>def dtdsNewton = { s, t, tR, k -> k * (tR - t) } assert dtdsNewton.maximumNumberOfParameters == 4

def dtds = dtdsNewton.rcurry(20, 0.07) assert dtds.maximumNumberOfParameters == 2

def tEulerH = eulerMapping.rcurry(dtds) { s, t -> s >= 100 } assert tEulerH.maximumNumberOfParameters == 3</lang>

Newton's Analytic Temperature Diffusion Solution (for comparison) <lang groovy>def tNewton = { s, s0, t0, tR, k ->

   tR + (t0 - tR) * Math.exp(k * (s0 - s))

} assert tNewton.maximumNumberOfParameters == 5

def tAnalytic = tNewton.rcurry(0, 100, 20, 0.07) assert tAnalytic.maximumNumberOfParameters == 1</lang>

Specific solutions for 3 step sizes (and initial time and temperature) <lang groovy>[10, 5, 2].each { h ->

   def tEuler = tEulerH.rcurry(h)
   assert tEuler.maximumNumberOfParameters == 2
   println """

STEP SIZE == ${h}

 time   analytic   euler   relative

(seconds) (°C) (°C) error

-------- -------- ---------"""

   tEuler(0, 100).each { BigDecimal s, tE -> 
       def tA = tAnalytic(s)
       def relError = ((tE - tA)/(tA - 20)).abs()
       printf('%5.0f    %8.4f %8.4f %9.6f\n', s, tA, tE, relError)
   }

}</lang>

Selected output

STEP SIZE == 10
  time   analytic   euler   relative
(seconds)  (°C)     (°C)     error
-------- -------- -------- ---------
    0    100.0000 100.0000  0.000000
   10     59.7268  44.0000  0.395874
   20     39.7278  27.2000  0.635032
   30     29.7965  22.1600  0.779513
   40     24.8648  20.6480  0.866798
   50     22.4158  20.1944  0.919529
   60     21.1996  20.0583  0.951386
   70     20.5957  20.0175  0.970631
   80     20.2958  20.0052  0.982257
   90     20.1469  20.0016  0.989281
  100     20.0730  20.0005  0.993524

STEP SIZE == 5
  time   analytic   euler   relative
(seconds)  (°C)     (°C)     error
-------- -------- -------- ---------
    0    100.0000 100.0000  0.000000
     ... yada, yada, yada ...
  100     20.0730  20.0145  0.801240

STEP SIZE == 2
  time   analytic   euler   relative
(seconds)  (°C)     (°C)     error
-------- -------- -------- ---------
    0    100.0000 100.0000  0.000000
     ... yada, yada, yada ...
  100     20.0730  20.0425  0.417918

Notice how the relative error in the Euler method sequences increases over time in spite of the fact that all three the Euler approximations and the analytic solution are approaching the same asymptotic limit of 20°C.

Notice also how smaller step size reduces the relative error in the approximation.

Haskell

<lang Haskell>import Text.Printf

type F = Double -> Double -> Double

euler :: F -> Double -> Double -> Double -> Double -> IO () euler f y0 a b h = do printf "%6.3f %6.3f\n" a y0 if a < b then euler f (y0 + (f a y0) * h) (a + h) b h else putStrLn "DONE"

newtonCooling :: F newtonCooling _ t = -0.07 * (t - 20)

main = euler newtonCooling 100 0 100 10 </lang> 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
100.000 20.000
DONE

Icon and Unicon

Translation of: Common Lisp

This solution works in both Icon and Unicon. It takes advantage of the proc procedure, which converts a string naming a procedure into a call to that procedure.

<lang Icon> invocable "newton_cooling" # needed to use the 'proc' procedure

procedure euler (f, y0, a, b, h)

 t := a
 y := y0
 until (t >= b) do {
   write (right(t, 4) || " " || left(y, 7))
   t +:= h
   y +:= h * (proc(f) (t, y)) # 'proc' applies procedure named in f to (t, y)
 }
 write ("DONE")

end

procedure newton_cooling (time, T)

 return -0.07 * (T - 20)

end

procedure main ()

 # generate data for all three step sizes [2, 5, 10]
 every (step_size := ![2,5,10]) do 
   euler ("newton_cooling", 100, 0, 100,  step_size)

end </lang>

Sample output:

J

Solution: <lang j>NB.*euler a Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and t=a..b and step size h. euler=: adverb define

'Y a b h'=. 4{. y
t=. a
while. b >: tincr=. {:t + h do.
  Y=. Y, (+ h * u) {:Y
  t=. t, tincr
end.
t,.Y

)

ncl=: _0.07 * -&20 NB. Newton's Cooling Law </lang> Example: <lang j> ncl euler 100 0 100 2 ... NB. output redacted for brevity

  ncl euler 100 0 100 5

... NB. output redacted for brevity

  ncl euler 100 0 100 10
 0     100
10      44
20    27.2
30   22.16
40  20.648
50 20.1944
60 20.0583
70 20.0175
80 20.0052
90 20.0016

100 20.0005</lang>

Java

<lang java> public class Euler {

 private static void euler (Callable f, double y0, int a, int b, int h) {
   int t = a;
   double y = y0;
   while (t < b) {
     System.out.println ("" + t + " " + y);
     t += h;
     y += h * f.compute (t, y);
   }
   System.out.println ("DONE");
 }

 public static void main (String[] args) {
   Callable cooling = new Cooling ();
   int[] steps = {2, 5, 10};
   for (int stepSize : steps) {
     System.out.println ("Step size: " + stepSize);
     euler (cooling, 100.0, 0, 100, stepSize);
   }
 }

}

// interface used so we can plug in alternative functions to Euler interface Callable {

 public double compute (int time, double t);

}

// class to implement the newton cooling equation class Cooling implements Callable {

 public double compute (int time, double t) {
   return -0.07 * (t - 20);
 }

} </lang>

Output for step = 10;

Step size: 10
0 100.0
10 43.99999999999999
20 27.199999999999996
30 22.159999999999997
40 20.648
50 20.194399999999998
60 20.05832
70 20.017496
80 20.0052488
90 20.00157464
DONE

Lua

<lang lua>T0 = 100 TR = 20 k = 0.07 delta_t = { 2, 5, 10 } n = 100

NewtonCooling = function( t ) return -k * ( t - TR ) end

function Euler( f, y0, n, h )

   local y = y0
   for x = 0, n, h do

print( "", x, y )

	y = y + h * f( y )
   end

end

for i = 1, #delta_t do

   print( "delta_t = ", delta_t[i] )
   Euler( NewtonCooling, T0, n, delta_t[i] )

end </lang>

Perl

<lang Perl>sub euler_method {

       my ($t0, $t1, $k, $step_size) = @_;
       my @results = ( [0, $t0] );

       for (my $s = $step_size; $s <= 100; $s += $step_size) {
               $t0 -= ($t0 - $t1) * $k * $step_size;
               push @results, [$s, $t0];
       }

       return @results;

}

sub analytical {

       my ($t0, $t1, $k, $time) = @_;
       return ($t0 - $t1) * exp(-$time * $k) + $t1

}

my ($T0, $T1, $k) = (100, 20, .07); my @r2 = grep { $_->[0] % 10 == 0 } euler_method($T0, $T1, $k, 2); my @r5 = grep { $_->[0] % 10 == 0 } euler_method($T0, $T1, $k, 5); my @r10 = grep { $_->[0] % 10 == 0 } euler_method($T0, $T1, $k, 10);

print "Time\t 2 err(%) 5 err(%) 10 err(%) Analytic\n", "-" x 76, "\n"; for (0 .. $#r2) {

       my $an = analytical($T0, $T1, $k, $r2[$_][0]);
       printf "%4d\t".("%9.3f" x 7)."\n",
               $r2 [$_][0],
               $r2 [$_][1], ($r2 [$_][1] / $an) * 100 - 100,
               $r5 [$_][1], ($r5 [$_][1] / $an) * 100 - 100,
               $r10[$_][1], ($r10[$_][1] / $an) * 100 - 100,
               $an;

} </lang>

Output:

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

Perl 6

<lang perl6>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 );
   }
   return @t_y;

}

my $COOLING_RATE = 0.07; my $AMBIENT_TEMP = 20; my $INITIAL_TEMP = 100; my $INITIAL_TIME = 0; my $FINAL_TIME = 100;

sub f ( $time, $temp ) {

   return -$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 * abs( $_ - $exact ) / $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%%');

}</lang>

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%

PicoLisp

<lang PicoLisp>(load "@lib/math.l")

(de euler (F Y A B H)

  (while (> B A)
     (prinl (round A) " " (round Y))
     (inc 'Y (*/ H (F A Y) 1.0))
     (inc 'A H) ) )

(de newtonCoolingLaw (A B)

  (*/ -0.07 (- B 20.) 1.0) )

(euler newtonCoolingLaw 100.0 0 100.0 2.0) (euler newtonCoolingLaw 100.0 0 100.0 5.0) (euler newtonCoolingLaw 100.0 0 100.0 10.0)</lang> 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.018
80.000 20.005
90.000 20.002

PureBasic

<lang PureBasic>Define.d Prototype.d Func(Time, t)

Procedure.d Euler(*F.Func, y0, a, b, h)

 Protected y=y0, t=a
 While t<=b
   PrintN(RSet(StrF(t,3),7)+" "+RSet(StrF(y,3),7))
   y + h * *F(t,y)
   t + h
 Wend

EndProcedure

Procedure.d newtonCoolingLaw(Time, t)

 ProcedureReturn -0.07*(t-20)

EndProcedure

If OpenConsole()

 Euler(@newtonCoolingLaw(), 100, 0, 100, 2)
 Euler(@newtonCoolingLaw(), 100, 0, 100, 5)
 Euler(@newtonCoolingLaw(), 100, 0, 100,10)
 
 Print(#CRLF$ + #CRLF$ + "Press ENTER to exit"): Input()
 CloseConsole()

EndIf</lang>

...
 85.000  20.053
 90.000  20.034
 95.000  20.022
100.000  20.014
  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
100.000  20.000

Python

Translation of: Common Lisp

<lang python>def euler(f,y0,a,b,h): t,y = a,y0 while t < b: print "%6.3f %6.3f" % (t,y) t += h y += h * f(t,y)

def newtoncooling(time, temp): return -0.07 * (temp - 20)

euler(newtoncooling,100,0,100,10) </lang> 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

Ruby

Translation of: Python

<lang ruby>def euler(y0,a,b,h, &block)

 t,y = a,y0
 while t < b
   puts "%6.3f %6.3f" % [t,y]
   t += h
   y += h * block.call(t,y)
 end

end

euler(100,0,100,10) {|time, temp| -0.07 * (temp - 20) } </lang> 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

Scala

<lang scala> object App{

 def main(args : Array[String]) = {
   
   def cooling( step : Int ) = {
     eulerStep( (step , y) => {-0.07 * (y - 20)} , 
       100.0,0,100,step)
   }
   cooling(10)
   cooling(5)
   cooling(2)
 }
 def eulerStep( func : (Int,Double) => Double,y0 : Double,
   begin : Int, end : Int , step : Int) = {
   
   println("Step size: %s".format(step))
   
   var current : Int = begin
   var y : Double = y0
   while( current <= end){
     println( "%d %.5f".format(current,y))
     current += step
     y += step * func(current,y)
   }
   
   println("DONE")
 }

} </lang>

Output for step = 10;

Step size: 10
0 100.00000
10 44.00000
20 27.20000
30 22.16000
40 20.64800
50 20.19440
60 20.05832
70 20.01750
80 20.00525
90 20.00157
DONE

Tcl

Translation of: C++

<lang tcl>proc euler {f y0 a b h} {

   puts "computing $f over \[$a..$b\], step $h"
   set y [expr {double($y0)}]
   for {set t [expr {double($a)}]} {$t < $b} {set t [expr {$t + $h}]} {

puts [format "%.3f\t%.3f" $t $y] set y [expr {$y + $h * double([$f $t $y])}]

   }
   puts "done"

}</lang> Demonstration with the Newton Cooling Law: <lang tcl>proc newtonCoolingLaw {time temp} {

   expr {-0.07 * ($temp - 20)}

}

euler newtonCoolingLaw 100 0 100 2 euler newtonCoolingLaw 100 0 100 5 euler newtonCoolingLaw 100 0 100 10</lang> End of output:

...
computing newtonCoolingLaw over [0..100], step 10
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