Euler method
From Rosetta Code
You are encouraged to solve this task according to the task description, using any language you may know.
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
- y(t0) = y0
To get a numeric solution, we replace the derivative on the LHS with a finite difference approximation:
then solve for y(t + h):
which is the same as
The iterative solution rule is then:
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(t0) = T0 cools down in an environment of temperature TR:
or
It says that the cooling rate 
of the object is proportional to the current temperature difference ΔT = (T(t) − TR) to the surrounding environment.
The analytical solution, which we will compare to the numerical approximation, is
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 T0 shall be 100 °C, the room temperature TR 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.
[edit] Ada
The solution is generic, usable for any floating point type. The package specification:
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;
The function Solve returns the solution of the differential equation for each of N+1 points, starting from the point T0. The implementation:
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;
The test program:
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;
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
[edit] ALGOL 68
Note: This specimen retains the original D coding style.#
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)
)
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
[edit] BBC BASIC
PROCeuler("-0.07*(y-20)", 100, 0, 100, 2)
PROCeuler("-0.07*(y-20)", 100, 0, 100, 5)
PROCeuler("-0.07*(y-20)", 100, 0, 100, 10)
END
DEF PROCeuler(df$, y, a, b, s)
LOCAL t, @%
@% = &2030A
t = a
WHILE t <= b
PRINT t, y
y += s * EVAL(df$)
t += s
ENDWHILE
ENDPROC
Output:
0.000 100.000
2.000 88.800
4.000 79.168
6.000 70.884
8.000 63.761
10.000 57.634
...
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
[edit] C
#include <stdio.h>output
#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;
}
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
[edit] C++
#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);
}
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
[edit] C#
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);
}
}
}
}
[edit] 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);
}
}
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
[edit] Common 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)
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
[edit] 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);
}
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
[edit] Euler Math Toolbox
>function dgleuler (f,x,y0) ...
$ y=zeros(size(x)); y[1]=y0;
$ for i=2 to cols(y);
$ y[i]=y[i-1]+f(x[i-1],y[i-1])*(x[i]-x[i-1]);
$ end;
$ return y;
$endfunction
>function f(x,y) := -k*(y-TR)
>k=0.07; TR=20; TS=100;
>x=0:1:100; dgleuler("f",x,TS)[-1]
20.0564137335
>x=0:2:100; dgleuler("f",x,TS)[-1]
20.0424631834
>TR+(TS-TR)*exp(-k*TS)
20.0729505572
>x=0:5:100; plot2d(x,dgleuler("f",x,TS)); ...
> plot2d(x,TR+(TS-TR)*exp(-k*x),>add,color=red);
>ode("f",x,TS)[-1] // Euler default solver LSODA
20.0729505568
>adaptiverunge("f",x,TS)[-1] // Adaptive Runge Method
20.0729505572
[edit] F#
let euler f (h : float) t0 y0 =
(t0, y0)
|> Seq.unfold (fun (t, y) -> Some((t,y), ((t + h), (y + h * (f t y)))))
let newtonCoolíng _ y = -0.07 * (y - 20.0)
[<EntryPoint>]
let main argv =
let f = newtonCoolíng
let a = 0.0
let y0 = 100.0
let b = 100.0
let h = 10.0
(euler newtonCoolíng h a y0)
|> Seq.takeWhile (fun (t,_) -> t <= b)
|> Seq.iter (printfn "%A")
0
Output for the above (step size 10)
(0.0, 100.0) (10.0, 44.0) (20.0, 27.2) (30.0, 22.16) (40.0, 20.648) (50.0, 20.1944) (60.0, 20.05832) (70.0, 20.017496) (80.0, 20.0052488) (90.0, 20.00157464) (100.0, 20.00047239)
[edit] 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
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
[edit] 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()
}
}
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
[edit] 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: xn and yn 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, x0 and y0 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.
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
Specific Euler Method Solution for the "Temperature Diffusion" Problem (with Newton's derivative formula and constants for environment temperature and object conductivity given)
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
Newton's Analytic Temperature Diffusion Solution (for comparison)
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
Specific solutions for 3 step sizes (and initial time and temperature)
[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)
}
}
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.
[edit] Haskell
import Text.Printf
euler :: (Num a, Ord a) => (a -> a -> a) -> a -> a -> a -> a -> [(a,a)]
euler f y0 a b h =
(a, y0) :
if a < b
then euler f (y0 + (f a y0) * h) (a + h) b h
else []
newtonCooling :: Double -> Double -> Double
newtonCooling _ t = -0.07 * (t - 20)
main = do
mapM_ (uncurry $ printf "%6.3f %6.3f\n") $ euler newtonCooling 100 0 100 10
putStrLn "DONE"
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
[edit] Icon and Unicon
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.
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
Sample output:
0 100 10 44.0 20 27.2 30 22.16 40 20.648 50 20.1944 60 20.0583 70 20.0174 80 20.0052 90 20.0015 DONE
[edit] J
Solution:
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
'Y0 a b h'=. 4{. y
t=. i.@>:&.(%&h) b - a
Y=. (+ h * u)^:(<#t) Y0
t,.Y
)
ncl=: _0.07 * -&20 NB. Newton's Cooling Law
Example:
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
[edit] 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);
}
}
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
[edit] 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
[edit] Mathematica
Better methods for differential equation solving are built into Mathematica, so the typical user would omit the Method and StartingStepSize options in the code below. However since the task requests Eulers method, here is the bad solution...
euler[step_, val_] := NDSolve[{T'[t] == -0.07 (T[t] - 20), T[0] == 100}, T, {t, 0, 100}, Method -> "ExplicitEuler", StartingStepSize -> step][[1, 1, 2]][val]
- Output:
euler[2, 100] 20.0425 euler[5, 100] 20.0145 euler[10, 100] 20.0005
[edit] МК-61/52
П2 С/П П3 С/П П4 ПП 19 ИП3 * ИП4
+ П4 С/П ИП2 ИП3 + П2 БП 05 ...
... ... ... ... ... ... ... ... ... В/О
Instead of dots typed calculation program equation f(u, t), where the arguments are t = Р2, u = Р4.
Input: Initial time С/П Time step С/П Initial value С/П.
The result is displayed on the indicator.
[edit] Objeck
class EulerMethod {
T0 : static : Float;
TR : static : Float;
k : static : Float;
delta_t : static : Float[];
n : static : Float;
function : Main(args : String[]) ~ Nil {
T0 := 100;
TR := 20;
k := 0.07;
delta_t := [2.0, 5.0, 10.0];
n := 100;
f := NewtonCooling(Float) ~ Float;
for(i := 0; i < delta_t->Size(); i+=1;) {
IO.Console->Print("delta_t = ")->PrintLine(delta_t[i]);
Euler(f, T0, n->As(Int), delta_t[i]);
};
}
function : native : NewtonCooling(t : Float) ~ Float {
return -1 * k * (t-TR);
}
function : native : Euler(f : (Float) ~ Float, y : Float, n : Int, h : Float) ~ Nil {
for(x := 0; x<=n; x+=h;) {
IO.Console->Print("\t")->Print(x)->Print("\t")->PrintLine(y);
y += h * f(y);
};
}
}
Output:
delta_t = 2
0 100
2 88.8
4 79.168
6 70.88448
...
delta_t = 10
0 100
10 44
20 27.2
30 22.16
40 20.648
[edit] Perl
sub euler_method {Output:
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;
}
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
[edit] Perl 6
sub euler ( &f, $y0, $a, $b, $h ) {Output:
my $y = $y0;
my @t_y;
for $a, * + $h ... * > $b -> $t {
@t_y[$t] = $y;
$y += $h * f( $t, $y );
}
return @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 ) {
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%%');
}
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%
[edit] 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)
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
[edit] PL/I
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;
Only the final two outputs are shown, for brevity.
H= 5; t By formula By Euler 0 100.0000000000 100.0000000000 5 76.3750457764 72.0000000000 10 59.7268257141 53.7999992371 15 47.9950218201 41.9700012207 20 39.7277565002 34.2805023193 25 33.9019165039 29.2823257446 30 29.7965145111 26.0335121155 35 26.9034862518 23.9217834473 40 24.8648052216 22.5491600037 45 23.4281692505 21.6569538116 50 22.4157905579 21.0770206451 55 21.7023792267 20.7000637054 60 21.1996459961 20.4550418854 65 20.8453769684 20.2957763672 70 20.5957260132 20.1922550201 75 20.4198017120 20.1249656677 80 20.2958297729 20.0812282562 85 20.2084674835 20.0527992249 90 20.1469039917 20.0343189240 95 20.1035213470 20.0223064423 100 20.0729503632 20.0144996643 H= 10; t By formula By Euler 0 100.0000000000 100.0000000000 10 59.7268257141 44.0000000000 20 39.7277565002 27.2000007629 30 29.7965145111 22.1599998474 40 24.8648052216 20.6480007172 50 22.4157905579 20.1944007874 60 21.1996459961 20.0583209991 70 20.5957260132 20.0174961090 80 20.2958297729 20.0052490234 90 20.1469039917 20.0015754700 100 20.0729503632 20.0004730225
[edit] 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
... 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
[edit] 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)
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
[edit] Racket
The ODE solver:
(define (ODE-solve f init
#:x-max x-max
#:step h
#:method (method euler))
(reverse
(iterate-while (λ (x . y) (<= x x-max)) (method f h) init)))
It uses the default integration method euler, defined separately.
(define (euler F h)
(λ (x y) (list (+ x h) (+ y (* h (F x y))))))
A general-purpose procedure which evalutes a given function f repeatedly starting with argument x, while all results satisfy a predicate test. Returns a list of iterations.
(define (iterate-while test f x)
(let next ([result x]
[list-of-results '()])
(if (apply test result)
(next (apply f result) (cons result list-of-results))
list-of-results)))
Textual output:
> (define (newton-cooling t T)
(* -0.07 (- T 20)))
> (ODE-solve newton-cooling '(0 100) #:x-max 100 #:step 10)
'((0 100)
(10 44.)
(20 27.2)
(30 22.16)
(40 20.648)
(50 20.1944)
(60 20.05832)
(70 20.017496)
(80 20.0052488)
(90 20.00157464)
(100 20.000472392))
Plotting results:
> (require plot)
> (plot
(map (λ (h c)
(lines
(ODE-solve newton-cooling '(0 100) #:x-max 100 #:step h)
#:color c #:label (format "h=~a" h)))
'(10 5 1)
'(red blue black))
#:legend-anchor 'top-right)
High modularity of the program allows to implement very different solution metods. For example 2-nd order Runge-Kutta method:
(define (RK2 F h)
(λ (x y)
(list (+ x h) (+ y (* h (F (+ x (* 1/2 h))
(+ y (* 1/2 h (F x y)))))))))
Two-step Adams–Bashforth method
(define (adams F h)
(case-lambda
; first step using Runge-Kutta method
[(x y) (append ((RK2 F h) x y) (list (F x y)))]
[(x y f′)
(let ([f (F x y)])
(list (+ x h) (+ y (* 3/2 h f) (* -1/2 h f′)) f))]))
Adaptive one-step method modifier using absolute accuracy ε
(define ((adaptive method ε) F h0)
(case-lambda
[(x y) (((adaptive method ε) F h0) x y h0)]
[(x y h)
(match-let* ([(list x0 y0) ((method F h) x y)]
[(list x1 y1) ((method F (/ h 2)) x y)]
[(list x1 y1) ((method F (/ h 2)) x1 y1)]
[τ (abs (- y1 y0))]
[h′ (if (< τ ε) (min h h0) (* 0.9 h (/ ε τ)))])
(list x1 (+ y1 τ) (* 2 h′)))]))
Comparison of different integration methods
> (define (solve-newton-cooling-by m)
(ODE-solve newton-cooling '(0 100)
#:x-max 100 #:step 10 #:method m))
> (plot
(list
(function (λ (t) (+ 20 (* 80 (exp (* -0.07 t))))) 0 100
#:color 'black #:label "analytical")
(lines (solve-newton-cooling-by euler)
#:color 'red #:label "Euler")
(lines (solve-newton-cooling-by RK2)
#:color 'blue #:label "Runge-Kutta")
(lines (solve-newton-cooling-by adams)
#:color 'purple #:label "Adams")
(points (solve-newton-cooling-by (adaptive euler 0.5))
#:color 'red #:label "Adaptive Euler")
(points (solve-newton-cooling-by (adaptive RK2 0.5))
#:color 'blue #:label "Adaptive Runge-Kutta"))
#:legend-anchor 'top-right)
See also Runge-Kutta method#Racket
[edit] 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) }
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
[edit] Run BASIC
x = euler(-0.07,-20, 100, 0, 100, 2)
x = euler-0.07,-20, 100, 0, 100, 5)
x = euler(-0.07,-20, 100, 0, 100, 10)
end
FUNCTION euler(da,db, y, a, b, s)
print "===== da:";da;" db:";db;" y:";y;" a:";a;" b:";b;" s:";s;" ==================="
t = a
WHILE t <= b
PRINT t;chr$(9);y
y = y + s * (da * (y + db))
t = t + s
WEND
END FUNCTION
===== da:-0.07 db:-20 y:100 a:0 b:100 s:2 =================== 0 100 2 88.8 4 79.168 6 70.88448 8 63.7606528 10 57.6341614 12 52.3653788 14 47.8342258 ...... ===== da:-0.07 db:-20 y:100 a:0 b:100 s:10 =================== 0 100 10 44.0 20 27.2 30 22.16 40 20.648 50 20.1944 60 20.05832 70 20.017496 80 20.0052488
[edit] 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")
}
}
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
[edit] Smalltalk
ODESolver>>eulerOf: f init: y0 from: a to: b step: h
| t y |
t := a.
y := y0.
[ t < b ]
whileTrue: [
Transcript
show: t asString, ' ' , (y printShowingDecimalPlaces: 3);
cr.
t := t + h.
y := y + (h * (f value: t value: y)) ]
ODESolver new eulerOf: [:time :temp| -0.07 * (temp - 20)] init: 100 from: 0 to: 100 step: 10
Transcript:
0 100.000 10 44.000 20 27.200 30 22.160 40 20.648 50 20.194 60 20.058 70 20.017 80 20.005 90 20.002
[edit] 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"
}
Demonstration with the Newton Cooling Law:
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
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
[edit] XPL0
include c:\cxpl\codes; \intrinsic 'code' declarations
proc Euler(Step); \Display cooling temperatures using Euler's method
int Step;
int Time; real Temp;
[Text(0, "Step "); IntOut(0, Step); Text(0, " ");
Time:= 0; Temp:= 100.0;
repeat if rem(Time/10) = 0 then RlOut(0, Temp);
Temp:= Temp + float(Step) * (-0.07*(Temp-20.0));
Time:= Time + Step;
until Time > 100;
CrLf(0);
];
real Time, Temp;
[Format(6,0); \display time heading
Text(0, "Time ");
Time:= 0.0;
while Time <= 100.1 do \(.1 avoids possible rounding error)
[RlOut(0, Time);
Time:= Time + 10.0;
];
CrLf(0);
Format(3,2); \display cooling temps using differential eqn.
Text(0, "Dif eq "); \ dTemp(time)/dtime = -k*Temp
Time:= 0.0;
while Time <= 100.1 do
[Temp:= 20.0 + (100.0-20.0) * Exp(-0.07*Time);
RlOut(0, Temp);
Time:= Time + 10.0;
];
CrLf(0);
Euler(2); \display cooling temps for various time steps
Euler(5);
Euler(10);
]
Output:
Time 0 10 20 30 40 50 60 70 80 90 100 Dif eq 100.00 59.73 39.73 29.80 24.86 22.42 21.20 20.60 20.30 20.15 20.07 Step 2 100.00 57.63 37.70 28.33 23.92 21.84 20.87 20.41 20.19 20.09 20.04 Step 5 100.00 53.80 34.28 26.03 22.55 21.08 20.46 20.19 20.08 20.03 20.01 Step 10 100.00 44.00 27.20 22.16 20.65 20.19 20.06 20.02 20.01 20.00 20.00
