Euler method: Difference between revisions

Content added Content deleted

Inline

Revision as of 22:10, 27 April 2011

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.

An 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++

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

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
DONE

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)(F f, double y0, double a, double b, double h) {

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

}

/// Example: Newton's cooling law auto newtonCoolingLaw(double time, double t) {

   return -0.07 * (t - 20);

}

void main() {

   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

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>

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>

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

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

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