Problem Statement

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: with an initial value To get a numeric solution, we replace the derivative on the   LHS   with a finite difference approximation: then solve for

y ( t + h )

{\displaystyle y(t+h)}

: which is the same as The iterative solution rule is then: 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}}

: or

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

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: and to compare with the analytical solution.

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

Let's start with the solution:

test: procedure options (main); /* 3 December 2012 */

   declare (x, y, z) float;
   declare (T0 initial (100), Tr initial (20)) float;
   declare k float initial (0.07);
   declare t fixed binary;
   declare h fixed binary;

   x, y, z = T0;
   /* Step size is 2 seconds */
   h = 2;
   put skip data (h);
   put skip list ('  t    By formula', 'By Euler');
   do t = 0 to 100 by 2;
      put skip edit (t, Tr + (T0 - Tr)/exp(k*t), x) (f(3), 2 f(17,10));
      x = x + h*f(t, x);
   end;

   /* Step size is 5 seconds */
   h = 5;
   put skip data (h);
   put skip list ('  t    By formula', 'By Euler');
   do t = 0 to 100 by 5;
      put skip edit ( t, Tr + (T0 - Tr)/exp(k*t), y) (f(3), 2 f(17,10));
      y = y + h*f(t, y);
   end;

   /* Step size is 10 seconds */
   h = 10;
   put skip data (h);
   put skip list ('  t    By formula', 'By Euler');
   do t = 0 to 100 by 10;
      put skip edit (t, Tr + (T0 - Tr)/exp(k*t), z) (f(3), 2 f(17,10));
      z = z + h*f(t, z);
   end;

f: procedure (dummy, T) returns (float);
   declare dummy  fixed binary;
   declare T float;

   return ( -k*(T - Tr) );
end f;

end test;