Euler method
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
To get a numeric solution, we replace the derivative on the LHS with a finite difference approximation:
then solve for :
which is the same as
The iterative solution rule is then:
where 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 cools down in an environment of temperature :
or
It says that the cooling rate of the object is proportional to the current temperature difference to the surrounding environment.
The analytical solution, which we will compare to the numerical approximation, is
- 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 shall be 100 °C
- room temperature shall be 20 °C
- cooling constant 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.
F euler(f, y0, a, b, h)
V t = a
V y = y0
L t <= b
print(‘#2.3 #2.3’.format(t, y))
t += h
y += h * f(t, y)
V newtoncooling = (time, temp) -> -0.07 * (temp - 20)
euler(newtoncooling, 100.0, 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.017 80.000 20.005 90.000 20.002 100.000 20.000
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
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
Which is
begin % Euler's method %
% Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and t=a..b and the step size h. %
real procedure euler ( real procedure f; real value y0, a, b, h ) ;
begin
real y, t;
y := y0;
t := a;
while t < b do begin
write( r_format := "A", r_w := 8, r_d := 4, s_w := 0, t, ": ", y );
y := y + ( h * f(t, y) );
t := t + h
end while_t_lt_b ;
write( "done" );
y
end euler ;
% Example: Newton's cooling law %
real procedure newtonCoolingLaw ( real value time, t ) ; -0.07 * (t - 20);
euler( newtonCoolingLaw, 100, 0, 100, 10 )
end.- Output:
0.0000: 100.0000 10.0000: 44.0000 20.0000: 27.2000 30.0000: 22.1600 40.0000: 20.6480 50.0000: 20.1944 60.0000: 20.0583 70.0000: 20.0175 80.0000: 20.0052 90.0000: 20.0015 done
euler: function [f, y0, a, b, h][
[t,y]: @[a, y0]
while [t < b][
print [to :string .format:".3f" t, to :string .format:".3f" y]
t: t + h
y: y + h * call f @[t,y]
]
]
newtoncooling: function [ti, te]->
(neg 0.07) * te - 20
euler 'newtoncooling 100.0 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.017 80.000 20.005 90.000 20.002
The following program's output should be fed to Gnuplot, which will produce a PNG (using either the font that is specified by the ATS program or a fallback substitute). You can run the program with a command such as this: patscc -g -O2 -std=gnu2x -DATS_MEMALLOC_LIBC euler_method_task.dats && ./a.out | gnuplot
All data requiring allocation is allocated as linear types, and so is not leaked. That is, roughly speaking: it requires no garbage collection.
#include "share/atspre_staload.hats"
staload UN = "prelude/SATS/unsafe.sats"
#define NIL list_vt_nil ()
#define :: list_vt_cons
(* Approximate y(t) in dy/dt=f(t,y), y(a)=y0, t going from a to b with
positive step size h. This implementation of euler_method requires
f to be a unboxed linear closure. *)
extern fn {tk : tkind}
euler_method (f : &(g0float tk, g0float tk) -<clo1> g0float tk,
y0 : g0float tk,
a : g0float tk,
b : g0float tk,
h : g0float tk) : List1_vt @(g0float tk, g0float tk)
implement {tk}
euler_method (f, y0, a, b, h) =
let
typedef point_pair = @(g0float tk, g0float tk)
fun
loop (f : &(g0float tk, g0float tk) -<clo1> g0float tk,
t : g0float tk,
y : g0float tk,
point_pairs : List0_vt point_pair)
: List1_vt point_pair =
let
val point_pairs = @(t, y) :: point_pairs
in
if b <= t then
reverse<point_pair> point_pairs
else
loop (f, t + h, y + (h * f (t, y)), point_pairs)
end
in
loop (f, a, y0, NIL)
end
fun {tk : tkind}
write_point_pairs
(outf : FILEref,
point_pairs : !List0_vt @(g0float tk, g0float tk))
: void =
case+ point_pairs of
| NIL => ()
| (@(t, y) :: tl) =>
begin
fprint_val<g0float tk> (outf, t);
fprint! (outf, " ");
fprint_val<g0float tk> (outf, y);
fprintln! (outf);
write_point_pairs (outf, tl)
end
implement
main0 () =
let
(* Implement f as a stack-allocated linear closure. *)
var f =
lam@ (t : double, y : double) : double => ~0.07 * (y - 20.0)
val data2 = euler_method<dblknd> (f, 100.0, 0.0, 100.0, 2.0)
and data5 = euler_method<dblknd> (f, 100.0, 0.0, 100.0, 5.0)
and data10 = euler_method<dblknd> (f, 100.0, 0.0, 100.0, 10.0)
val outf = stdout_ref
in
fprintln! (outf, "set encoding utf8");
fprintln! (outf, "set term png size 1000,750 font 'RTF Amethyst Pro,16'");
fprintln! (outf, "set output 'newton-cooling-ATS.png'");
fprintln! (outf, "set grid");
fprintln! (outf, "set title 'Newton’s Law of Cooling'");
fprintln! (outf, "set xlabel 'Elapsed time (seconds)'");
fprintln! (outf, "set ylabel 'Temperature (Celsius)'");
fprintln! (outf, "set xrange [0:100]");
fprintln! (outf, "set yrange [15:100]");
fprintln! (outf, "y(x) = 20.0 + (80.0 * exp (-0.07 * x))");
fprintln! (outf, "plot y(x) with lines title 'Analytic solution', \\");
fprintln! (outf, " '-' with linespoints title 'Euler method, step size 2s', \\");
fprintln! (outf, " '-' with linespoints title 'Step size 5s', \\");
fprintln! (outf, " '-' with linespoints title 'Step size 10s'");
write_point_pairs (outf, data2);
fprintln! (outf, "e");
write_point_pairs (outf, data5);
fprintln! (outf, "e");
write_point_pairs (outf, data10);
fprintln! (outf, "e");
free data2;
free data5;
free data10
end- Output:
100 HOME
110 PRINT "Time ";
120 FOR S = 0 TO 100.1 STEP 10
130 PRINT S; " ";
140 NEXT S
150 PRINT
160 PRINT "Dif eq ";
170 FOR S = 0 TO 100.1 STEP 10
180 LET T = 20+(100-20)*EXP(-.07*S)
190 PRINT LEFT$(STR$(T),5); " ";
200 NEXT S
210 PRINT
220 LET P = 2 : GOSUB 260
230 LET P = 5 : GOSUB 260
240 LET P = 10 : GOSUB 260
250 END
260 REM Euler(paso)
270 LET S = 0
280 LET T = 100
290 PRINT "Step ";P; " ";
300 FOR S = 0 TO 100 STEP P
310 IF S - (S/10) * 10 = 0 THEN PRINT LEFT$(STR$(T),5); " ";
320 LET T = T+(P)*(-.07*(T-20))
330 IF S > 100 THEN GOTO 350
340 NEXT S
350 PRINT
360 RETURN
print "Time ";
tiempo = 0.0
while tiempo <= 100.1
print rjust(string(tiempo), 5); " ";
tiempo += 10.0
end while
print
print "Dif eq ";
tiempo = 0.0
while tiempo <= 100.1
temperatura = 20.0 + (100.0-20.0) * exp(-0.07*tiempo)
print rjust(left(string(temperatura), 5), 5); " ";
tiempo += 10.0
end while
print
call Euler(2)
call Euler(5)
call Euler(10)
end
subroutine Euler(paso)
tiempo = 0
temperatura = 100.0
print "Step "; rjust(string(paso), 2); " ";
while tiempo <= 100
if (tiempo mod 10) = 0 then print rjust(left(string(temperatura), 5), 5); " ";
temperatura += float(paso) * (-0.07*(temperatura-20.0))
tiempo += paso
end while
print
end subroutine
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
Same code as GW-BASIC
precision 4
let s = 2
gosub euler
let s = 5
gosub euler
let s = 10
gosub euler
end
sub euler
cls
cursor 1, 1
wait
print "step: ", s
let b = 100
let y = 100
for t = 0 to b step s
print t, " : ", y
let y = y + s * (-0.07 * (y - 20))
gosub delay
next t
alert "step ", s, " finished"
return
sub delay
let w = clock
do
wait
loop clock < w + 200
return
- Output:
step: 2100 88.8000 79.1680 70.8845 63.7607 57.6342 52.3654 47.8342 43.9374 40.5862 37.7041 35.2255 33.0939 31.2608 29.6843 28.3285 27.1625 26.1597 25.2973 24.5557 23.9179 23.3694 22.8977 22.4920 22.1431 21.8431 21.5851 21.3632 21.1724 21.0083 20.8671 20.7457 20.6413 20.5515 20.4743 20.4079 20.3508 20.3017 20.2595 20.2232 20.1920 20.1651 20.1420 20.1221 20.1050 20.0903 20.0777 20.0668 20.0574 20.0494 20.0425
step: 5 100 72 53.8000 41.9700 34.2805 29.2823 26.0335 23.9218 22.5492 21.6570 21.0771 20.7001 20.4551 20.2958 20.1923 20.1250 20.0813 20.0528 20.0343 20.0223 20.0145
step: 10
100 44 27.2000 22.1600 20.6480 20.1944 20.0583 20.0175 20.0053 20.0016 20.0005
'Freebasic .9
'Custom rounding
#define round(x,N) Rtrim(Rtrim(Left(Str((x)+(.5*Sgn((x)))/(10^(N))),Instr(Str((x)+(.5*Sgn((x)))/(10^(N))),".")+(N)),"0"),".")
#macro Euler(fn,_y,min,max,h,printoption)
Print "Step ";#h;":":Print
Print "time","Euler"," Analytic"
If printoption<>"print" Then Print "Data omitted ..."
Scope
Dim As Double temp=(min),y=(_y)
Do
If printoption="print" Then Print temp,round(y,3),20+80*Exp(-0.07*temp)
y=y+(h)*(fn)
temp=temp+(h)
Loop Until temp>(max)
Print"________________"
Print
End Scope
#endmacro
Euler(-.07*(y-20),100,0,100,2,"don't print")
Euler(-.07*(y-20),100,0,100,5,"print")
Euler(-.07*(y-20),100,0,100,10,"print")
Sleep
outputs (steps 5 and 10)
Step 2: time Euler Analytic Data omitted ... ________________ Step 5: time Euler Analytic 0 100 100 5 72 76.37504717749707 10 53.8 59.72682430331276 15 41.97 47.99501992889243 20 34.281 39.72775711532852 25 29.282 33.90191547603561 30 26.034 29.79651426023855 35 23.922 26.90348691994964 40 22.549 24.86480501001743 45 21.657 23.42817014936322 50 21.077 22.41579067378548 55 20.7 21.70237891507017 60 20.455 21.19964614563822 65 20.296 20.84537635070821 70 20.192 20.59572664567395 75 20.125 20.41980147193451 80 20.081 20.29582909731863 85 20.053 20.20846724147268 90 20.034 20.14690438216231 95 20.022 20.10352176843727 100 20.014 20.07295055724436 ________________ Step 10: time Euler Analytic 0 100 100 10 44 59.72682430331276 20 27.2 39.72775711532852 30 22.16 29.79651426023855 40 20.648 24.86480501001743 50 20.194 22.41579067378548 60 20.058 21.19964614563822 70 20.017 20.59572664567395 80 20.005 20.29582909731863 90 20.002 20.14690438216231 100 20 20.07295055724436 ________________
Public Sub Main()
Dim tiempo As Float, temperatura As Float
Print "Time ";
tiempo = 0.0
While tiempo <= 100.1
Print Format$(tiempo, "#######");
tiempo += 10.0
Wend
Print
Print "Dif eq ";
tiempo = 0.0
While tiempo <= 100.1
temperatura = 20.0 + (100.0 - 20.0) * Exp(-0.07 * tiempo)
Print Format$(temperatura, "####.#0");
tiempo += 10.0
Wend
Print
Euler(2)
Euler(5)
Euler(10)
End
Public Sub Euler(paso As Integer)
Dim tiempo As Integer = 0
Dim temperatura As Float = 100.0
Print "Step "; Format$(paso, "##"); " ";
Do While tiempo <= 100
If (tiempo Mod 10) = 0 Then Print Format$(temperatura, "####.#0");
temperatura += (paso) * (-0.07 * (temperatura - 20.0))
tiempo += paso
Loop
Print
End Sub
100 CLS
110 PRINT "Time ";
120 FOR TIEMPO = 0 TO 100.1 STEP 10
130 PRINT USING " ###";TIEMPO;
140 NEXT TIEMPO
150 PRINT
160 PRINT "Dif eq ";
170 FOR TIEMPO = 0 TO 100.1 STEP 10
180 TEMPERATURA = 20+(100-20)*EXP(-.07*TIEMPO)
190 PRINT USING "###.##";TEMPERATURA;
200 NEXT TIEMPO
210 PRINT
220 PASO = 2 : GOSUB 260
230 PASO = 5 : GOSUB 260
240 PASO = 10 : GOSUB 260
250 END
260 REM Euler(paso)
270 TIEMPO = 0
280 TEMPERATURA = 100
290 PRINT USING "Step ## ";PASO;
300 FOR TIEMPO = 0 TO 100 STEP PASO
310 IF (TIEMPO MOD 10) = 0 THEN PRINT USING "###.##";TEMPERATURA;
320 TEMPERATURA = TEMPERATURA+(PASO)*(-.07*(TEMPERATURA-20))
330 IF TIEMPO > 100 THEN EXIT FOR
340 NEXT TIEMPO
350 PRINT
360 RETURN
Same code as GW-BASIC
DECLARE SUB Euler (paso AS INTEGER)
CLS
PRINT "Time ";
tiempo = 0!
WHILE tiempo <= 100.1
PRINT USING "######"; tiempo;
tiempo = tiempo + 10!
WEND
PRINT
PRINT "Dif eq ";
tiempo = 0!
WHILE tiempo <= 100.1
temperatura = 20! + (100! - 20!) * EXP(-.07 * tiempo)
PRINT USING "###.##"; temperatura;
tiempo = tiempo + 10!
WEND
PRINT
Euler (2)
Euler (5)
Euler (10)
END
SUB Euler (paso AS INTEGER)
tiempo = 0
temperatura = 100!
PRINT USING "Step ## "; paso;
DO WHILE tiempo <= 100
IF (tiempo MOD 10) = 0 THEN PRINT USING "###.##"; temperatura;
temperatura = temperatura + paso * (-.07 * (temperatura - 20!))
tiempo = tiempo + paso
LOOP
PRINT
END SUB
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
SUB euler (paso)
LET tiempo = 0
LET temperatura = 100
PRINT USING "Step ## ": paso;
DO WHILE tiempo <= 100
IF (REMAINDER(tiempo,10)) = 0 THEN PRINT USING "###.##": temperatura;
LET temperatura = temperatura+paso*(-.07*(temperatura-20))
LET tiempo = tiempo+paso
LOOP
PRINT
END SUB
PRINT "Time ";
LET tiempo = 0
DO WHILE tiempo <= 100.1
PRINT USING "######": tiempo;
LET tiempo = tiempo+10
LOOP
PRINT
PRINT "Dif eq ";
LET tiempo = 0
DO WHILE tiempo <= 100.1
LET temperatura = 20+(100-20)*EXP(-.07*tiempo)
PRINT USING "###.##": temperatura;
LET tiempo = tiempo+10
LOOP
PRINT
CALL Euler (2)
CALL Euler (5)
CALL Euler (10)
END
- Output:
Same as QBasic entry.
PROGRAM "Euclidean rhythm"
VERSION "0.0001"
IMPORT "xma"
DECLARE FUNCTION Entry ()
DECLARE FUNCTION Euler (paso)
FUNCTION Entry ()
PRINT "Time ";
tiempo! = 0.0
DO WHILE tiempo! <= 100.1
PRINT FORMAT$ ("#######", tiempo!);
tiempo! = tiempo! + 10.0
LOOP
PRINT
PRINT "Dif eq ";
tiempo! = 0.0
DO WHILE tiempo! <= 100.1
temperatura! = 20.0 + (100.0 - 20.0) * EXP(-0.07 * tiempo!)
PRINT FORMAT$ ("####.##", temperatura!);
tiempo! = tiempo! + 10.0
LOOP
PRINT
Euler(2)
Euler(5)
Euler(10)
END FUNCTION
FUNCTION Euler (paso)
tiempo! = 0
temperatura! = 100.0
PRINT FORMAT$ ("Step ## ", paso);
DO WHILE tiempo! <= 100
IF (tiempo! MOD 10) = 0 THEN PRINT FORMAT$ ("####.##", temperatura!);
temperatura! = temperatura! + SINGLE(paso) * (-0.07 * (temperatura! - 20.0))
tiempo! = tiempo! + paso
LOOP
PRINT
END FUNCTION
END PROGRAM
print "Time ";
tiempo = 0.0
while tiempo <= 100.1
print tiempo using "#######";
tiempo = tiempo + 10.0
wend
print
print "Dif eq ";
tiempo = 0.0
while tiempo <= 100.1
temperatura = 20.0 + (100.0-20.0) * exp(-0.07*tiempo)
print temperatura using "####.##";
tiempo = tiempo + 10.0
wend
print
Euler_(2)
Euler_(5)
Euler_(10)
end
sub Euler_(paso)
local tiempo, temperatura
tiempo = 0
temperatura = 100.0
print "Step ", paso using "##", " ";
while tiempo <= 100
if mod(tiempo, 10) = 0 print temperatura using "####.##";
temperatura = temperatura + (paso) * (-0.07*(temperatura-20.0))
tiempo = tiempo + paso
end while
print
end sub
#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;
}
output
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
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);
}
}
}
}
#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
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
(ns newton-cooling
(:gen-class))
(defn euler [f y0 a b h]
"Euler's Method.
Approximates y(time) in y'(time)=f(time,y) with y(a)=y0 and t=a..b and the step size h."
(loop [t a
y y0
result []]
(if (<= t b)
(recur (+ t h) (+ y (* (f (+ t h) y) h)) (conj result [(double t) (double y)]))
result)))
(defn newton-coolling [t temp]
"Newton's cooling law, f(t,T) = -0.07*(T-20)"
(* -0.07 (- temp 20)))
; Run for case h = 10
(println "Example output")
(doseq [q (euler newton-coolling 100 0 100 10)]
(println (apply format "%.3f %.3f" q)))
- Output:
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 100.000 20.000
The following is in the Managed COBOL dialect:
DELEGATE-ID func.
PROCEDURE DIVISION USING VALUE t AS FLOAT-LONG
RETURNING ret AS FLOAT-LONG.
END DELEGATE.
CLASS-ID. MainClass.
78 T0 VALUE 100.0.
78 TR VALUE 20.0.
78 k VALUE 0.07.
01 delta-t INITIALIZE ONLY STATIC
FLOAT-LONG OCCURS 3 VALUES 2.0, 5.0, 10.0.
78 n VALUE 100.
METHOD-ID NewtonCooling STATIC.
PROCEDURE DIVISION USING VALUE t AS FLOAT-LONG
RETURNING ret AS FLOAT-LONG.
COMPUTE ret = - k * (t - TR)
END METHOD.
METHOD-ID Main STATIC.
DECLARE f AS TYPE func
SET f TO METHOD self::NewtonCooling
DECLARE delta-t-len AS BINARY-LONG
MOVE delta-t::Length TO delta-t-len
PERFORM VARYING i AS BINARY-LONG FROM 1 BY 1
UNTIL i > delta-t-len
DECLARE elt AS FLOAT-LONG = delta-t (i)
INVOKE TYPE Console::WriteLine("delta-t = {0:F4}", elt)
INVOKE self::Euler(f, T0, n, elt)
END-PERFORM
END METHOD.
METHOD-ID Euler STATIC.
PROCEDURE DIVISION USING VALUE f AS TYPE func, y AS FLOAT-LONG,
n AS BINARY-LONG, h AS FLOAT-LONG.
PERFORM VARYING x AS BINARY-LONG FROM 0 BY h UNTIL x >= n
INVOKE TYPE Console::WriteLine("x = {0:F4}, y = {1:F4}", x, y)
COMPUTE y = y + h * RUN f(y)
END-PERFORM
END METHOD.
END CLASS.
Example output:
delta-t = 10.0000 x = 0.0000, y = 100.0000 x = 10.0000, y = 44.0000 x = 20.0000, y = 27.2000 x = 30.0000, y = 22.1600 x = 40.0000, y = 20.6480 x = 50.0000, y = 20.1944 x = 60.0000, y = 20.0583 x = 70.0000, y = 20.0175 x = 80.0000, y = 20.0052 x = 90.0000, y = 20.0016
;; '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 (+ t h))
(y y0 (+ 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)
;; slightly more idiomatic Common Lisp version
(defun newton-cooling (time temperature)
"Newton's cooling law, f(t,T) = -0.07*(T-20)"
(declare (ignore time))
(* -0.07 (- temperature 20)))
(defun euler (f y0 a b h)
"Euler's Method.
Approximates y(time) in y'(time)=f(time,y) with y(a)=y0 and t=a..b and the step size h."
(loop for time from a below b by h
for y = y0 then (+ y (* h (funcall f time y)))
do (format t "~6,3F ~6,3F~%" time y)))
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
import std.stdio, std.range, std.traits;
/// 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) @safe
if (isCallable!F && __traits(compiles, { real r = f(0.0, 0.0); })) {
double y = y0;
foreach (immutable t; iota(a, b, h)) {
writefln("%.3f %.3f", t, y);
y += h * f(t, y);
}
"done".writeln;
}
void main() {
/// Example: Newton's cooling law.
enum newtonCoolingLaw = (in double time, in double t)
pure nothrow @safe @nogc => -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
import 'dart:math';
import "dart:io";
const double k = 0.07;
const double initialTemp = 100.0;
const double finalTemp = 20.0;
const int startTime = 0;
const int endTime = 100;
void ivpEuler(double Function(double, double) function, double initialValue, int step) {
stdout.write(' Step ${step.toString().padLeft(2)}: ');
var y = initialValue;
for (int t = startTime; t <= endTime; t += step) {
if (t % 10 == 0) {
stdout.write(y.toStringAsFixed(3).padLeft(7));
}
y += step * function(t.toDouble(), y);
}
print('');
}
void analytic() {
stdout.write(' Time: ');
for (int t = startTime; t <= endTime; t += 10) {
stdout.write(t.toString().padLeft(7));
}
stdout.write('\nAnalytic: ');
for (int t = startTime; t <= endTime; t += 10) {
var temp = finalTemp + (initialTemp - finalTemp) * exp(-k * t);
stdout.write(temp.toStringAsFixed(3).padLeft(7));
}
print('');
}
double cooling(double t, double temp) {
return -k * (temp - finalTemp);
}
void main() {
analytic();
ivpEuler(cooling, initialTemp, 2);
ivpEuler(cooling, initialTemp, 5);
ivpEuler(cooling, initialTemp, 10);
}
- Output:
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
DuckDB-defined functions are not first-class citizens, and a function must have been declared before it is referenced, but it is permissible to redefine a function, so in the following we illustrate that a dummy declaration for cooling/2 can be provided separately, thus affording a degree of modularity.
# Declaration for the sake of euler_table()
create or replace function cooling(x,y) as NULL;
# Euler method assuming cooling/2 is available
create or replace function euler_table(x_start, y_start, x2, h) as table (
with recursive cte as (
select x_start::DOUBLE as x, y_start::DOUBLE as y
union all
select x + h as x,
y + h*cooling(x, y)
from cte
where (x < x2 and x_start < x2) or
(x > x2 and x_start > x2)
)
from cte
);
create or replace function euler(x_start, y_start, x2, h) as (
select last(y order by x) from euler_method(x_start, y_start, x2, h)
);
## Example:
create or replace function cooling(x,y) as
-0.07 * (y - 20);
select h, euler(0,100, 100, h) from unnest([1,2,5,10,20]) _(h);
- Output:
┌───────┬───────────────────────┐ │ h │ euler(0, 100, 100, h) │ │ int32 │ double │ ├───────┼───────────────────────┤ │ 1 │ 20.05641373347389 │ │ 2 │ 20.0424631833732 │ │ 5 │ 20.01449963666907 │ │ 10 │ 20.000472392 │ │ 20 │ 19.180799999999998 │ └───────┴───────────────────────┘
TR = 20
K = -0.07
ytxt = 95
drawgrid
glinewidth 0.3
gtextsize 3
#
func newton_cooling temp .
return K * (temp - TR)
.
proc draw_euler y t t2 step col .
gcolor col
gtext 80 ytxt "step: " & step
ytxt -= 5
gpenup
while t <= t2
glineto t y
t += step
y += step * newton_cooling y
.
.
func analytic T0 t .
return TR + (T0 - TR) * pow 2.71828 (K * t)
.
proc draw_analytic a b .
gcolor 009
gtext 80 ytxt "analytic"
ytxt -= 5
gpenup
for t = a to b
y = analytic 100 t
glineto t y
.
.
draw_euler 100 0 100 10 900
draw_euler 100 0 100 5 555
draw_euler 100 0 100 2 090
draw_analytic 0 100defmodule Euler do
def method(_, _, t, b, _) when t>b, do: :ok
def method(f, y, t, b, h) do
:io.format "~7.3f ~7.3f~n", [t,y]
method(f, y + h * f.(t,y), t + h, b, h)
end
end
f = fn _time, temp -> -0.07 * (temp - 20) end
Enum.each([10, 5, 2], fn step ->
IO.puts "\nStep = #{step}"
Euler.method(f, 100.0, 0.0, 100.0, step)
end)
- Output:
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 100.000 20.000 Step = 5 0.000 100.000 5.000 72.000 10.000 53.800 15.000 41.970 20.000 34.280 25.000 29.282 30.000 26.034 35.000 23.922 40.000 22.549 45.000 21.657 50.000 21.077 55.000 20.700 60.000 20.455 65.000 20.296 70.000 20.192 75.000 20.125 80.000 20.081 85.000 20.053 90.000 20.034 95.000 20.022 100.000 20.014 Step = 2 0.000 100.000 2.000 88.800 4.000 79.168 6.000 70.884 8.000 63.761 10.000 57.634 12.000 52.365 14.000 47.834 16.000 43.937 18.000 40.586 20.000 37.704 22.000 35.226 24.000 33.094 26.000 31.261 28.000 29.684 30.000 28.328 32.000 27.163 34.000 26.160 36.000 25.297 38.000 24.556 40.000 23.918 42.000 23.369 44.000 22.898 46.000 22.492 48.000 22.143 50.000 21.843 52.000 21.585 54.000 21.363 56.000 21.172 58.000 21.008 60.000 20.867 62.000 20.746 64.000 20.641 66.000 20.551 68.000 20.474 70.000 20.408 72.000 20.351 74.000 20.302 76.000 20.259 78.000 20.223 80.000 20.192 82.000 20.165 84.000 20.142 86.000 20.122 88.000 20.105 90.000 20.090 92.000 20.078 94.000 20.067 96.000 20.057 98.000 20.049 100.000 20.042
-module(euler).
-export([main/0, euler/5]).
cooling(_Time, Temperature) ->
(-0.07)*(Temperature-20).
euler(_, Y, T, _, End) when End == T ->
io:fwrite("\n"),
Y;
euler(Func, Y, T, Step, End) ->
if
T rem 10 == 0 ->
io:fwrite("~.3f ",[float(Y)]);
true ->
ok
end,
euler(Func, Y + Step * Func(T, Y), T + Step, Step, End).
analytic(T, End) when T == End ->
io:fwrite("\n"),
T;
analytic(T, End) ->
Y = (20 + 80 * math:exp(-0.07 * T)),
io:fwrite("~.3f ", [Y]),
analytic(T+10, End).
main() ->
io:fwrite("Analytic:\n"),
analytic(0, 100),
io:fwrite("Step 2:\n"),
euler(fun cooling/2, 100, 0, 2, 100),
io:fwrite("Step 5:\n"),
euler(fun cooling/2, 100, 0, 5, 100),
io:fwrite("Step 10:\n"),
euler(fun cooling/2, 100, 0, 10, 100),
ok.
- Output:
Analytic: 100.000 59.727 39.728 29.797 24.865 22.416 21.200 20.596 20.296 20.147 Step 2: 100.000 57.634 37.704 28.328 23.918 21.843 20.867 20.408 20.192 20.090 Step 5: 100.000 53.800 34.280 26.034 22.549 21.077 20.455 20.192 20.081 20.034 Step 10: 100.000 44.000 27.200 22.160 20.648 20.194 20.058 20.017 20.005 20.002 ok
>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.0729505572let 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)
USING: formatting fry io kernel locals math math.ranges
sequences ;
IN: rosetta-code.euler-method
:: euler ( quot y! a b h -- )
a b h <range> [
:> t
t y "%7.3f %7.3f\n" printf
t y quot call h * y + y!
] each ; inline
: cooling ( t y -- x ) nip 20 - -0.07 * ;
: euler-method-demo ( -- )
2 5 10 [ '[ [ cooling ] 100 0 100 _ euler ] call nl ] tri@ ;
MAIN: euler-method-demo
- 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
: newton-cooling-law ( f: temp -- f: temp' )
20e f- -0.07e f* ;
: euler ( f: y0 xt step end -- )
1+ 0 do
cr i . fdup f.
fdup over execute
dup s>f f* f+
dup +loop
2drop fdrop ;
100e ' newton-cooling-law 2 100 euler cr
100e ' newton-cooling-law 5 100 euler cr
100e ' newton-cooling-law 10 100 euler cr
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
Specialised to the cooling function. We produce an array of the temperature at each step subtracted from the analytically determined temperature (so we are computing the error).
let analytic(t0: f64) (time: f64): f64 =
20.0 + (t0 - 20.0) * f64.exp(-0.07*time)
let cooling(_time: f64) (temperature: f64): f64 =
-0.07 * (temperature-20.0)
let main(t0: f64) (a: f64) (b: f64) (h: f64): []f64 =
let steps = i32.f64 ((b-a)/h)
let temps = replicate steps 0.0
let (_,temps) = loop (t,temps)=(t0,temps) for i < steps do
let x = a + f64.i32 i * h
let temps[i] = f64.abs(t-analytic t0 x)
in (t + h * cooling x t,
temps)
in temps
FB 7.0.34 macOS 14.7.2 Sonoma
include "NSLog.incl"
double local fn Cooling( t as double, temp as double )
return -0.07 * (temp - 20 )
end fn = 0.0
void local fn Euler( y as double, stepSize as int, endTime as double )
int t = 0
NSLog( @"\n Step %2d: \b", stepSize )
while ( t <= endTime )
if ( t % 10 == 0 ) then NSLog( @" %7.3f \b", y )
y += stepSize * fn Cooling( t, y )
t += stepSize
wend
end fn
void local fn Analytic
NSLog( @" Time: \b" )
for int t1 = 0 to 100 step 10
NSLog( @" %8d\b", t1 )
next
NSLog( @"\nAnalytic: \b" )
for int t2 = 0 to 100 step 10
NSLog( @"%9.3f\b", 20 + 80 * exp(-0.07 * t2 ) )
next
end fn
fn Analytic
fn Euler( 100.0, 2.0, 100.0 )
fn Euler( 100.0, 5.0, 100.0 )
fn Euler( 100.0, 10.0, 100.0 )
HandleEvents- Output:
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
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
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.
Modular solution which separates the solver and a method. Moreover it works on a given mesh which can be irregular.
-- the solver
dsolveBy _ _ [] _ = error "empty solution interval"
dsolveBy method f mesh x0 = zip mesh results
where results = scanl (method f) x0 intervals
intervals = zip mesh (tail mesh)
It is better to use strict Data.List.scanl' in the solver but avoiding highlighting problems we leave lazy scanl function.
Some possible methods:
-- 1-st order Euler
euler f x (t1,t2) = x + (t2 - t1) * f t1 x
-- 2-nd order Runge-Kutta
rk2 f x (t1,t2) = x + h * f (t1 + h/2) (x + h/2*f t1 x)
where h = t2 - t1
-- 4-th order Runge-Kutta
rk4 f x (t1,t2) = x + h/6 * (k1 + 2*k2 + 2*k3 + k4)
where k1 = f t1 x
k2 = f (t1 + h/2) (x + h/2*k1)
k3 = f (t1 + h/2) (x + h/2*k2)
k4 = f (t1 + h) (x + h*k3)
h = t2 - t1
Graphical output, using EasyPlot:
import Graphics.EasyPlot
newton t temp = -0.07 * (temp - 20)
exactSolution t = 80*exp(-0.07*t)+20
test1 = plot (PNG "euler1.png")
[ Data2D [Title "Step 10", Style Lines] [] sol1
, Data2D [Title "Step 5", Style Lines] [] sol2
, Data2D [Title "Step 1", Style Lines] [] sol3
, Function2D [Title "exact solution"] [Range 0 100] exactSolution ]
where sol1 = dsolveBy euler newton [0,10..100] 100
sol2 = dsolveBy euler newton [0,5..100] 100
sol3 = dsolveBy euler newton [0,1..100] 100
test2 = plot (PNG "euler2.png")
[ Data2D [Title "Euler"] [] sol1
, Data2D [Title "RK2"] [] sol2
, Data2D [Title "RK4"] [] sol3
, Function2D [Title "exact solution"] [Range 0 100] exactSolution ]
where sol1 = dsolveBy euler newton [0,10..100] 100
sol2 = dsolveBy rk2 newton [0,10..100] 100
sol3 = dsolveBy rk4 newton [0,10..100] 100
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.
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
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
For NGN/K
// make global variables with :: because while is only locally scoped
tn1::100
tn::100
// to store solution
sol::()
// these are constant, so only local scope
k:-0.07
time:100
steps:10
h:time%steps
tr:20
// applying x+h to 0 outside to end while loop
// bind to a variable to surpress output
var: {x<100}{tn1::tn+h*k*(tn-tr);sol::sol,tn;tn::tn1;x+h}\0
sol
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
// Function that takes differential-equation, initial condition,
// ending x, and step size as parameters
function eulersMethod(f, x1, y1, x2, h) {
// Header
console.log("\tX\t|\tY\t");
console.log("------------------------------------");
// Initial Variables
var x=x1, y=y1;
// While we're not done yet
// Both sides of the OR let you do Euler's Method backwards
while ((x<x2 && x1<x2) || (x>x2 && x1>x2)) {
// Print what we have
console.log("\t" + x + "\t|\t" + y);
// Calculate the next values
y += h*f(x, y)
x += h;
}
return y;
}
function cooling(x, y) {
return -0.07 * (y-20);
}
eulersMethod(cooling, 0, 100, 100, 10);
# euler_method takes a filter (df), initial condition
# (x1,y1), ending x (x2), and step size as parameters;
# it emits the y values at each iteration.
# df must take [x,y] as its input.
def euler_method(df; x1; y1; x2; h):
h as $h
| [x1, y1]
| recurse( if ((.[0] < x2 and x1 < x2) or
(.[0] > x2 and x1 > x2)) then
[ (.[0] + $h), (.[1] + $h*df) ]
else empty
end )
| .[1] ;
# We could now solve the task by writing for each step-size, $h
# euler_method(-0.07 * (.[1]-20); 0; 100; 100; $h)
# but for clarity, we shall define a function named "cooling":
# [x,y] is input
def cooling: -0.07 * (.[1]-20);
# The following solves the task:
# (2,5,10) | [., [ euler_method(cooling; 0; 100; 100; .) ] ]For brevity, we modify euler_method so that it only shows the final value of y:
def euler_solution(df; x1; y1; x2; h):
def recursion(exp): reduce recurse(exp) as $x (.; $x);
h as $h
| [x1, y1]
| recursion( if ((.[0] < x2 and x1 < x2) or
(.[0] > x2 and x1 > x2)) then
[ (.[0] + $h), (.[1] + $h*df) ]
else empty
end )
| .[1] ;Example:
(1,2,5,10,20) | [., [ euler_solution(cooling; 0; 100; 100; .) ] ]- Output:
$ jq -M -n -c -f Euler_method.jq
[1,[20.05641373347389]]
[2,[20.0424631833732]]
[5,[20.01449963666907]]
[10,[20.000472392]]
[20,[19.180799999999998]]
euler(f::Function, T::Number, t0::Int, t1::Int, h::Int) = collect(begin T += h * f(T); T end for t in t0:h:t1)
# Prints a series of arbitrary values in a tabular form, left aligned in cells with a given width
tabular(width, cells...) = println(join(map(s -> rpad(s, width), cells)))
# prints the table according to the task description for h=5 and 10 sec
for h in (5, 10)
print("Step $h:\n\n")
tabular(15, "Time", "Euler", "Analytic")
t = 0
for T in euler(y -> -0.07 * (y - 20.0), 100.0, 0, 100, h)
tabular(15, t, round(T,digits=6), round(20.0 + 80.0 * exp(-0.07t), digits=6))
t += h
end
println()
end
- Output:
Step 5: Time Euler Analytic 0 72.0 100.0 5 53.8 76.375047 10 41.97 59.726824 15 34.2805 47.99502 20 29.282325 39.727757 25 26.033511 33.901915 30 23.921782 29.796514 35 22.549159 26.903487 40 21.656953 24.864805 45 21.077019 23.42817 50 20.700063 22.415791 55 20.455041 21.702379 60 20.295776 21.199646 65 20.192255 20.845376 70 20.124966 20.595727 75 20.081228 20.419801 80 20.052798 20.295829 85 20.034319 20.208467 90 20.022307 20.146904 95 20.0145 20.103522 100 20.009425 20.072951 Step 10: Time Euler Analytic 0 44.0 100.0 10 27.2 59.726824 20 22.16 39.727757 30 20.648 29.796514 40 20.1944 24.864805 50 20.05832 22.415791 60 20.017496 21.199646 70 20.005249 20.595727 80 20.001575 20.295829 90 20.000472 20.146904 100 20.000142 20.072951
// version 1.1.2
typealias Deriv = (Double) -> Double // only one parameter needed here
const val FMT = " %7.3f"
fun euler(f: Deriv, y: Double, step: Int, end: Int) {
var yy = y
print(" Step %2d: ".format(step))
for (t in 0..end step step) {
if (t % 10 == 0) print(FMT.format(yy))
yy += step * f(yy)
}
println()
}
fun analytic() {
print(" Time: ")
for (t in 0..100 step 10) print(" %7d".format(t))
print("\nAnalytic: ")
for (t in 0..100 step 10)
print(FMT.format(20.0 + 80.0 * Math.exp(-0.07 * t)))
println()
}
fun cooling(temp: Double) = -0.07 * (temp - 20.0)
fun main(args: Array<String>) {
analytic()
for (i in listOf(2, 5, 10))
euler(::cooling, 100.0, i, 100)
}
- Output:
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
Translation of Python
{def eulersMethod
{def eulersMethod.r
{lambda {:f :b :h :t :y}
{if {<= :t :b}
then {tr {td :t} {td {/ {round {* :y 1000}} 1000}}}
{eulersMethod.r :f :b :h {+ :t :h} {+ :y {* :h {:f :t :y}}}}
else}}}
{lambda {:f :y0 :a :b :h}
{table {eulersMethod.r :f :b :h :a :y0}}}}
{def cooling
{lambda {:time :temp}
{* -0.07 {- :temp 20}}}}
{eulersMethod cooling 100 0 100 10}
->
0 100
10 44
20 27.2
30 22.16
40 20.648
50 20.194
60 20.058
70 20.017
80 20.005
90 20.002
100 20
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
Build-in function with Euler:
with(Student[NumericalAnalysis]);
k := 0.07:
TR := 20:
Euler(diff(T(t), t) = -k*(T(t) - TR), T(0) = 100, t = 100, numsteps = 50); # step size = 2
Euler(diff(T(t), t) = -k*(T(t) - TR), T(0) = 100, t = 100, numsteps = 20); # step size = 5
Euler(diff(T(t), t) = -k*(T(t) - TR), T(0) = 100, t = 100, numsteps = 10); # step size = 10
- Output:
20.04 20.01 20.00
Hard-coded procedure:
f := y -> (-0.07) * (y - 20):
EulerMethod := proc(f, start_time, end_time, y0, h) # y0: initial value #h: step size
local cur, y, rate:
cur := start_time;
y := y0;
while cur <= end_time do
printf("%g %g\n", cur, y);
cur := cur + h;
rate := f(y);
y := y + h * rate;
end do;
return y;
end proc:
# step size = 2
printf("Step Size = %a\n", 2);
EulerMethod(f, 0, 100, 100, 2);
# step size = 5
printf("\nStep Size = %a\n", 5);
EulerMethod(f, 0, 100, 100, 5);
# step size = 10
printf("\nStep Size = %a\n", 10);
EulerMethod(f, 0, 100, 100, 10);
- Output:
Step Size = 2 0 100 2 88.8 4 79.168 6 70.8845 8 63.7607 10 57.6342 12 52.3654 14 47.8342 16 43.9374 18 40.5862 20 37.7041 22 35.2255 24 33.094 26 31.2608 28 29.6843 30 28.3285 32 27.1625 34 26.1598 36 25.2974 38 24.5558 40 23.918 42 23.3694 44 22.8977 46 22.492 48 22.1432 50 21.8431 52 21.5851 54 21.3632 56 21.1723 58 21.0082 60 20.867 62 20.7457 64 20.6413 66 20.5515 68 20.4743 70 20.4079 72 20.3508 74 20.3017 76 20.2594 78 20.2231 80 20.1919 82 20.165 84 20.1419 86 20.122 88 20.105 90 20.0903 92 20.0776 94 20.0668 96 20.0574 98 20.0494 100 20.0425 Step Size = 5 0 100 5 72 10 53.8 15 41.97 20 34.2805 25 29.2823 30 26.0335 35 23.9218 40 22.5492 45 21.657 50 21.077 55 20.7001 60 20.455 65 20.2958 70 20.1923 75 20.125 80 20.0812 85 20.0528 90 20.0343 95 20.0223 100 20.0145 Step Size = 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
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
clear all;close all;clc;
format longG;
% Main Script
for h = [5, 10]
fprintf('Step %d:\n\n', h);
tabular(15, 'Time', 'Euler', 'Analytic');
T0 = 100.0;
t0 = 0;
t1 = 100;
T = euler(@(T) -0.07 * (T - 20.0), T0, t0, t1, h);
for i = 1:length(T)
t = (i-1) * h;
analytic = 20.0 + 80.0 * exp(-0.07 * t);
tabular(15, t, round(T(i), 6), round(analytic, 6));
end
fprintf('\n');
end
function T = euler(f, T0, t0, t1, h)
% EULER A simple implementation of Euler's method for solving ODEs
% f - function handle for the derivative
% T0 - initial temperature
% t0, t1 - start and end times
% h - step size
T = T0;
for t = t0:h:t1
T(end+1) = T(end) + h * f(T(end));
end
end
function tabular(width, varargin)
% TABULAR Prints a series of values in a tabular form
% width - cell width
% varargin - variable number of arguments representing cells
for i = 1:length(varargin)
fprintf('%-*s', width, num2str(varargin{i}));
end
fprintf('\n');
end- Output:
Step 5: Time Euler Analytic 0 100 100 5 72 76.375 10 53.8 59.7268 15 41.97 47.995 20 34.2805 39.7278 25 29.2823 33.9019 30 26.0335 29.7965 35 23.9218 26.9035 40 22.5492 24.8648 45 21.657 23.4282 50 21.077 22.4158 55 20.7001 21.7024 60 20.455 21.1996 65 20.2958 20.8454 70 20.1923 20.5957 75 20.125 20.4198 80 20.0812 20.2958 85 20.0528 20.2085 90 20.0343 20.1469 95 20.0223 20.1035 100 20.0145 20.073 105 20.0094 20.0514 Step 10: Time Euler Analytic 0 100 100 10 44 59.7268 20 27.2 39.7278 30 22.16 29.7965 40 20.648 24.8648 50 20.1944 22.4158 60 20.0583 21.1996 70 20.0175 20.5957 80 20.0052 20.2958 90 20.0016 20.1469 100 20.0005 20.073 110 20.0001 20.0362
euler_method(f, y0, a, b, h):= block(
[t: a, y: y0, tg: [a], yg: [y0]],
unless t>=b do (
t: t + h,
y: y + f(t, y)*h,
tg: endcons(t, tg),
yg: endcons(y, yg)
),
[tg, yg]
);
/* initial temperature */
T0: 100;
/* environment of temperature */
Tr: 20;
/* the cooling constant */
k: 0.07;
/* end of integration */
tmax: 100;
/* analytical solution */
Tref(t):= Tr + (T0 - Tr)*exp(-k*t);
/* cooling rate */
dT(t, T):= -k*(T-Tr);
/* get numerical solution */
h: 10;
[tg, yg]: euler_method('dT, T0, 0, tmax, h);
/* plot analytical and numerical solution */
plot2d([Tref, [discrete, tg, yg]], ['t, 0, tmax],
[legend, "analytical", concat("h = ", h)],
[xlabel, "t / seconds"],
[ylabel, "Temperature / C"]);
П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.
import strutils
proc euler(f: proc (x,y: float): float; y0, a, b, h: float) =
var (t,y) = (a,y0)
while t < b:
echo formatFloat(t, ffDecimal, 3), " ", formatFloat(y, ffDecimal, 3)
t += h
y += h * f(t,y)
proc newtoncooling(time, temp: float): float =
-0.07 * (temp - 20)
euler(newtoncooling, 100.0, 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.017 80.000 20.005 90.000 20.002
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
Output is a PNG produced by Gnuplot, which is run as a child process. The program demonstrates both UTF-8 support and ipl.functional.
Note the string written with u" because it contains a non-ASCII character (the apostrophe).
import
io(open),
ipl.functional(_a, lambda)
$encoding UTF-8
procedure main ()
local f, plot, ty
local data2, data5, data10
# Newton's cooling law, f(t,Temp) = -0.07*(Temp-20)
f := lambda { -0.07 * (_a[2] - 20.0) }
data2 := euler_method (f, 100, 0, 100, 2)
data5 := euler_method (f, 100, 0, 100, 5)
data10 := euler_method (f, 100, 0, 100, 10)
plot := open ("gnuplot", "pw")
plot.write ("set encoding utf8")
plot.write ("set term png size 1000,750 font 'Fanwood Text,18'")
plot.write ("set output 'newton-cooling-OI.png'")
plot.write ("set grid")
plot.write (u"set title 'Newton’s Law of Cooling'")
plot.write ("set xlabel 'Elapsed time (seconds)'")
plot.write ("set ylabel 'Temperature (Celsius)'")
plot.write ("set xrange [0:100]")
plot.write ("set yrange [15:100]")
plot.write ("y(x) = 20.0 + (80.0 * exp (-0.07 * x))")
plot.write ("plot y(x) with lines title 'Analytic solution', \\")
plot.write (" '-' with linespoints title 'Euler method, step size 2s', \\")
plot.write (" '-' with linespoints title 'Step size 5s', \\")
plot.write (" '-' with linespoints title 'Step size 10s'")
every plot.write (ty := !data2 & ty[1] || " " || ty[2])
plot.write ("e")
every plot.write (ty := !data5 & ty[1] || " " || ty[2])
plot.write ("e")
every plot.write (ty := !data10 & ty[1] || " " || ty[2])
plot.write ("e")
plot.close()
end
# Approximate y(t) in dy/dt=f(t,y), y(a)=y0, t going from a to b with
# positive step size h.
procedure euler_method (f, y0, a, b, h)
local t, y, results
t := a
y := y0
results := [[t, y]]
while t + h <= b do
{
y +:= h * f(t, y)
t +:= h
put (results, [t, y])
}
return results
end- Output:
(* Euler integration by recurrence relation.
* Given a function, and stepsize, provides a function of (t,y) which
* returns the next step: (t',y'). *)
let euler f ~step (t,y) = ( t+.step, y +. step *. f t y )
(* newton_cooling doesn't use time parameter, so _ is a placeholder *)
let newton_cooling ~k ~tr _ y = -.k *. (y -. tr)
(* analytic solution for Newton cooling *)
let analytic_solution ~k ~tr ~t0 t = tr +. (t0 -. tr) *. exp (-.k *. t)
Using the above functions to produce the task results:
(* Wrapping up the parameters in a "cool" function: *)
let cool = euler (newton_cooling ~k:0.07 ~tr:20.)
(* Similarly for the analytic solution: *)
let analytic = analytic_solution ~k:0.07 ~tr:20. ~t0:100.
(* (Just a loop) Apply recurrence function on state, until some condition *)
let recur ~until f state =
let rec loop s =
if until s then ()
else loop (f s)
in loop state
(* 'results' generates the specified output starting from initial values t=0, temp=100C; ending at t=100s *)
let results fn =
Printf.printf "\t time\t euler\tanalytic\n%!";
let until (t,y) =
Printf.printf "\t%7.3f\t%7.3f\t%9.5f\n%!" t y (analytic t);
t >= 100.
in recur ~until fn (0.,100.)
results (cool ~step:10.)
results (cool ~step:5.)
results (cool ~step:2.)
Example output:
# results (cool ~step:10.);;
time euler analytic
0.000 100.000 100.00000
10.000 44.000 59.72682
20.000 27.200 39.72776
30.000 22.160 29.79651
40.000 20.648 24.86481
50.000 20.194 22.41579
60.000 20.058 21.19965
70.000 20.017 20.59573
80.000 20.005 20.29583
90.000 20.002 20.14690
100.000 20.000 20.07295
- : unit = ()
package main
import "core:fmt"
import "core:math"
ivp_euler ::proc(y:f32,step:i32,end:i32){
y:=y
t:=i32(0)
for ;t<=end;t+=step{
if(t%10==0){fmt.printfln("%.3f", y)}
y+=f32(step)*cooling(f32(t),y)
}
}
cooling :: proc(t:f32,temp:f32)-> f32{
return -0.07* (temp-20.0)
}
analytic :: proc(){
for t := 0; t <= 100; t += 10{
fmt.printfln("%.3f",20+80*math.exp_f32(f32(-0.07)*f32(t)))
}
}
main :: proc() {
fmt.println("Step: 2 Seconds")
ivp_euler(100.0, 2,100)
fmt.println("Step: 5 Seconds")
ivp_euler(100.0, 5, 100)
fmt.println("Step: 10 Seconds")
ivp_euler(100.0, 10, 100)
fmt.println("Analytic")
analytic()
}
- Output:
Step: 2 Seconds 100.000 57.634 37.704 28.328 23.918 21.843 20.867 20.408 20.192 20.090 20.042 Step: 5 Seconds 100.000 53.800 34.281 26.034 22.549 21.077 20.455 20.192 20.081 20.034 20.014 Step: 10 Seconds 100.000 44.000 27.200 22.160 20.648 20.194 20.058 20.017 20.005 20.002 20.000 Analytic 100.000 59.727 39.728 29.797 24.865 22.416 21.200 20.596 20.296 20.147 20.073
: euler(f, y, a, b, h)
| t |
a b h step: t [
System.Out t <<wjp(6, JUSTIFY_RIGHT, 3) " : " << y << cr
t y f perform h * y + ->y
] ;Usage :
: newtonCoolingLaw(t, y)
y 20 - -0.07 * ;
: test
euler(#newtonCoolingLaw, 100.0, 0.0, 100.0, 2)
euler(#newtonCoolingLaw, 100.0, 0.0, 100.0, 5)
euler(#newtonCoolingLaw, 100.0, 0.0, 100.0, 10) ;- Output:
....
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
See also Scheme.
The output is meant to be fed into Gnuplot. You can run the program like this: ol name-you-saved-the-program-as.scm | gnuplot
(Gnuplot may substitute a different font.)
(define (euler-method f y0 a b h)
;; Approximate y(t) in dy/dt=f(t,y), y(a)=y0, t going from a to b
;; with positive step size h. Produce a list of point pairs as
;; output.
(let loop ((t a)
(y y0)
(point-pairs '()))
(let ((point-pairs (cons (cons t y) point-pairs)))
(if (<= b t)
(reverse point-pairs)
(loop (+ t h) (+ y (* h (f t y))) point-pairs)))))
(define (newton-cooling-step t Temperature)
;; Newton's cooling law, with temperature in Celsius:
;;
;; f(t, Temperature) = -0.07*(Temperature - 20)
;;
(* -0.07 (- Temperature 20)))
(define data-for-stepsize=2
(euler-method newton-cooling-step 100.0 0.0 100.0 2.0))
(define data-for-stepsize=5
(euler-method newton-cooling-step 100.0 0.0 100.0 5.0))
(define data-for-stepsize=10
(euler-method newton-cooling-step 100.0 0.0 100.0 10.0))
(define (display-point-pairs point-pairs)
(let loop ((p point-pairs))
(if (pair? p)
(begin
(display (inexact (caar p)))
(display " ")
(display (inexact (cdar p)))
(newline)
(loop (cdr p))))))
(display "set encoding utf8") (newline)
(display "set term png size 1000,750 font 'Farao Book,16'") (newline)
(display "set output 'newton-cooling-Scheme.png'") (newline)
(display "set grid") (newline)
(display "set title 'Newton’s Law of Cooling'") (newline)
(display "set xlabel 'Elapsed time (seconds)'") (newline)
(display "set ylabel 'Temperature (Celsius)'") (newline)
(display "set xrange [0:100]") (newline)
(display "set yrange [15:100]") (newline)
(display "y(x) = 20.0 + (80.0 * exp (-0.07 * x))") (newline)
(display "plot y(x) with lines title 'Analytic solution', \\") (newline)
(display " '-' with linespoints title 'Euler method, step size 2s', \\") (newline)
(display " '-' with linespoints title 'Step size 5s', \\") (newline)
(display " '-' with linespoints title 'Step size 10s'") (newline)
(display-point-pairs data-for-stepsize=2)
(display "e") (newline)
(display-point-pairs data-for-stepsize=5)
(display "e") (newline)
(display-point-pairs data-for-stepsize=10)
(display "e") (newline)
- Output:
Euler code for Free Pascal - Delphi mode. Apart from the function-pointer calling convention for the NewtonCooling method, this example is ISO-7185 standard Pascal.
{$mode delphi}
PROGRAM Euler;
TYPE TNewtonCooling = FUNCTION (t: REAL) : REAL;
CONST T0 : REAL = 100.0;
CONST TR : REAL = 20.0;
CONST k : REAL = 0.07;
CONST time : INTEGER = 100;
CONST step : INTEGER = 10;
CONST dt : ARRAY[0..3] of REAL = (1.0,2.0,5.0,10.0);
VAR i : INTEGER;
FUNCTION NewtonCooling(t: REAL) : REAL;
BEGIN
NewtonCooling := -k * (t-TR);
END;
PROCEDURE Euler(F: TNewtonCooling; y, h : REAL; n: INTEGER);
VAR i: INTEGER = 0;
BEGIN
WRITE('dt=',trunc(h):2,':');
REPEAT
IF (i mod 10 = 0) THEN WRITE(' ',y:2:3);
INC(i,trunc(h));
y := y + h * F(y);
UNTIL (i >= n);
WRITELN;
END;
PROCEDURE Sigma;
VAR t: INTEGER = 0;
BEGIN
WRITE('Sigma:');
REPEAT
WRITE(' ',(20 + 80 * exp(-0.07 * t)):2:3);
INC(t,step);
UNTIL (t>=time);
WRITELN;
END;
BEGIN
WRITELN('Newton cooling function: Analytic solution (Sigma) with 3 Euler approximations.');
WRITELN('Time: ',0:7,10:7,20:7,30:7,40:7,50:7,60:7,70:7,80:7,90:7);
Sigma;
FOR i := 1 to 3 DO
Euler(NewtonCooling,T0,dt[i],time);
END.
Output:
Newton cooling function: Analytic solution (Sigma) with 3 Euler approximations. Time: 0 10 20 30 40 50 60 70 80 90 Sigma: 100.000 59.727 39.728 29.797 24.865 22.416 21.200 20.596 20.296 20.147 dt= 2: 100.000 57.634 37.704 28.328 23.918 21.843 20.867 20.408 20.192 20.090 dt= 5: 100.000 53.800 34.280 26.034 22.549 21.077 20.455 20.192 20.081 20.034 dt=10: 100.000 44.000 27.200 22.160 20.648 20.194 20.058 20.017 20.005 20.002
procedure euler(f: (real) -> real; y0, a, b, h: real);
begin
var (t, y) := (a, y0);
write(' Step', h:3, ':');
while t <= b do
begin
if integer(t) mod 10 = 0 then write(y:8:3);
t += h;
y += h * f(y);
end;
writeln;
end;
procedure analytic();
begin
write(' Time:');
for var t := 0 to 100 step 10 do write(t:8);
writeln;
write('Analytic:');
for var t := 0 to 100 step 10 do
write(20.0 + 80.0 * exp(-0.07 * t):8:3);
println();
end;
function newtoncooling(temp: real): real := -0.07 * (temp - 20);
begin
analytic;
foreach var i in |2, 5, 10| do
euler(newtoncooling, 100.0, 0.0, 100.0, i)
end.
- Output:
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.281 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
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;
}
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
You can run this online here.
-- -- demo\rosetta\Euler_method.exw -- ============================= -- with javascript_semantics function ivp_euler(atom y, integer f, step, end_t) sequence res = {} for t=0 to end_t by step do if remainder(t,10)==0 then res &= y end if y += step * call_func(f,{t, y}) end for return res end function function analytic() sequence res = {} for t=0 to 100 by 10 do res &= 20 + 80 * exp(-0.07 * t) end for return res end function function cooling(atom /*t*/, temp) return -0.07 * (temp - 20) end function constant x = tagset(100,0,10), a = analytic(), e2 = ivp_euler(100,cooling,2,100), e5 = ivp_euler(100,cooling,5,100), e10 = ivp_euler(100,cooling,10,100) printf(1," Time: %s\n",{join(x,fmt:="%7d")}) printf(1,"Analytic: %s\n",{join(a,fmt:="%7.3f")}) printf(1," Step 2: %s\n",{join(e2,fmt:="%7.3f")}) printf(1," Step 5: %s\n",{join(e5,fmt:="%7.3f")}) printf(1," Step 10: %s\n",{join(e10,fmt:="%7.3f")}) -- and a simple plot include pGUI.e include IupGraph.e function get_data(Ihandle /*graph*/) return {{"NAMES",{"analytical","h=2","h=5","h=10"}}, {x,a,CD_BLUE},{x,e2,CD_GREEN},{x,e5,CD_BLACK},{x,e10,CD_RED}} end function IupOpen() Ihandle graph = IupGraph(get_data,`RASTERSIZE=340x240,GRID=NO`) IupSetAttributes(graph,`XTICK=20,XMIN=0,XMAX=100,XMARGIN=25`) IupSetAttributes(graph,`YTICK=20,YMIN=20,YMAX=100`) IupShow(IupDialog(graph,`TITLE="Euler Method",MINSIZE=260x200`)) if platform()!=JS then IupMainLoop() IupClose() end if
- Output:
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
(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
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
local fmt = require "fmt"
local function euler(f, y, step, finish)
fmt.write(" Step %2d: ", step)
for t = 0, finish, step do
if t % 10 == 0 then fmt.write(" %7.3f", y) end
y += step * f(y)
end
print()
end
local function analytic()
io.write(" Time: ")
for t = 0, 100, 10 do fmt.write(" %7d", t) end
io.write("\nAnalytic: ")
for t = 0, 100, 10 do
fmt.write(" %7.3f", 20 + 80 * math.exp(-0.07 * t))
end
print()
end
local cooling = |temp| -> -0.07 * (temp - 20)
analytic()
for {2, 5, 10} as i do euler(cooling, 100, i, 100) end
- Output:
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
function euler (${f}, ${y}, $y0, $t0, $tEnd) {
function f-euler ($tn, $yn, $h) {
$yn + $h*(f $tn $yn)
}
function time ($t0, $h, $tEnd) {
$end = [MATH]::Floor(($tEnd - $t0)/$h)
foreach ($_ in 0..$end) { $_*$h + $t0 }
}
$time = time $t0 10 $tEnd
$time5 = time $t0 5 $tEnd
$time2 = time $t0 2 $tEnd
$yn10 = $yn5 = $yn2 = $y0
$i2 = $i5 = 0
foreach ($tn10 in $time) {
while($time2[$i2] -ne $tn10) {
$i2++
$yn2 = (f-euler $time2[$i2] $yn2 2)
}
while($time5[$i5] -ne $tn10) {
$i5++
$yn5 = (f-euler $time5[$i5] $yn5 5)
}
[pscustomobject]@{
t = "$tn10"
Analytical = "$("{0:N5}" -f (y $tn10))"
"Euler h = 2" = "$("{0:N5}" -f $yn2)"
"Euler h = 5" = "$("{0:N5}" -f $yn5)"
"Euler h = 10" = "$("{0:N5}" -f $yn10)"
"Error h = 2" = "$("{0:N5}" -f [MATH]::abs($yn2 - (y $tn10)))"
"Error h = 5" = "$("{0:N5}" -f [MATH]::abs($yn5 - (y $tn10)))"
"Error h = 10" = "$("{0:N5}" -f [MATH]::abs($yn10 - (y $tn10)))"
}
$yn10 = (f-euler $tn10 $yn10 10)
}
}
$k, $yr, $y0, $t0, $tEnd = 0.07, 20, 100, 0, 100
function f ($t, $y) {
-$k *($y - $yr)
}
function y ($t) {
$yr + ($y0 - $yr)*[MATH]::Exp(-$k*$t)
}
euler f y $y0 $t0 $tEnd | Format-Table -AutoSize
Output:
t Analytical Euler h = 2 Euler h = 5 Euler h = 10 Error h = 2 Error h = 5 Error h = 10 - ---------- ----------- ----------- ------------ ----------- ----------- ------------ 0 100.00000 100.00000 100.00000 100.00000 0.00000 0.00000 0.00000 10 59.72682 57.63416 53.80000 44.00000 2.09266 5.92682 15.72682 20 39.72776 37.70413 34.28050 27.20000 2.02363 5.44726 12.52776 30 29.79651 28.32850 26.03351 22.16000 1.46801 3.76300 7.63651 40 24.86481 23.91795 22.54916 20.64800 0.94685 2.31565 4.21681 50 22.41579 21.84311 21.07702 20.19440 0.57268 1.33877 2.22139 60 21.19965 20.86705 20.45504 20.05832 0.33260 0.74461 1.14133 70 20.59573 20.40788 20.19225 20.01750 0.18784 0.40347 0.57823 80 20.29583 20.19188 20.08123 20.00525 0.10395 0.21460 0.29058 90 20.14690 20.09027 20.03432 20.00157 0.05664 0.11259 0.14533 100 20.07295 20.04246 20.01450 20.00047 0.03049 0.05845 0.07248
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
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 100.000 20.000
t:0
time:100
k:-0.07
tr:20
sol:()
steps:10
h:time%steps
tn:100
while [t<time;tn1:tn + h*(k)*(tn-tr);sol:sol,tn1;tn:tn1;t:t+h]
sol
sol 44 27.2 22.16 20.648 20.1944 20.05832 20.0175 20.00525 20.00157 20.00047
euler <- function(f, y0, a, b, h)
{
t <- a
y <- y0
while (t < b)
{
cat(sprintf("%6.3f %6.3f\n", t, y))
t <- t + h
y <- y + h*f(t, y)
}
}
newtoncooling <- function(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
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
(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%
Rebol [
title: "Rosetta code: Euler method"
file: %Euler_method.r3
url: https://rosettacode.org/wiki/Euler_method
]
euler: func[
"Euler's method for Newton's cooling law, with formatted comparison"
step [integer!] "time step in seconds"
precision [decimal!] "rounding granularity for printed numbers"
][
print ["^/STEP:" step]
print "Time Euler Analytic"
print "-------------------------"
;; Initialize:
;; b: upper time bound (seconds) — here we reuse the initial temperature value 100
;; y: the current Euler-approximated temperature; initial condition T(0) = 100
b: y: 100
for time 0 b step [
printf [-3 " | " 9 "| "] reduce [
time
round/to y precision
round/to (20 + (80 * exp (-0.07 * time))) precision
]
;; Euler step update:
;; This applies the discrete forward Euler update using the ODE's RHS
y: y + (step * (-0.07 * (y - 20)))
]
]
;; Run three experiments with different step sizes
euler 2 0.0001
euler 5 0.0001
euler 10 0.0001
- Output:
STEP: 2 Time Euler Analytic ------------------------- 0 | 100.0 | 100.0 2 | 88.8 | 89.5487 4 | 79.168 | 80.4627 6 | 70.8845 | 72.5637 8 | 63.7607 | 65.6967 10 | 57.6342 | 59.7268 12 | 52.3654 | 54.5368 14 | 47.8342 | 50.0249 16 | 43.9374 | 46.1024 18 | 40.5862 | 42.6923 20 | 37.7041 | 39.7278 22 | 35.2255 | 37.1505 24 | 33.094 | 34.9099 26 | 31.2608 | 32.9621 28 | 29.6843 | 31.2687 30 | 28.3285 | 29.7965 32 | 27.1625 | 28.5167 34 | 26.1598 | 27.404 36 | 25.2974 | 26.4368 38 | 24.5558 | 25.5959 40 | 23.918 | 24.8648 42 | 23.3694 | 24.2293 44 | 22.8977 | 23.6767 46 | 22.492 | 23.1964 48 | 22.1432 | 22.7788 50 | 21.8431 | 22.4158 52 | 21.5851 | 22.1002 54 | 21.3632 | 21.8258 56 | 21.1723 | 21.5873 58 | 21.0082 | 21.3799 60 | 20.867 | 21.1996 62 | 20.7457 | 21.0429 64 | 20.6413 | 20.9067 66 | 20.5515 | 20.7882 68 | 20.4743 | 20.6852 70 | 20.4079 | 20.5957 72 | 20.3508 | 20.5179 74 | 20.3017 | 20.4502 76 | 20.2594 | 20.3914 78 | 20.2231 | 20.3403 80 | 20.1919 | 20.2958 82 | 20.165 | 20.2572 84 | 20.1419 | 20.2236 86 | 20.122 | 20.1944 88 | 20.105 | 20.169 90 | 20.0903 | 20.1469 92 | 20.0776 | 20.1277 94 | 20.0668 | 20.111 96 | 20.0574 | 20.0965 98 | 20.0494 | 20.0839 100 | 20.0425 | 20.073 STEP: 5 Time Euler Analytic ------------------------- 0 | 100.0 | 100.0 5 | 72.0 | 76.375 10 | 53.8 | 59.7268 15 | 41.97 | 47.995 20 | 34.2805 | 39.7278 25 | 29.2823 | 33.9019 30 | 26.0335 | 29.7965 35 | 23.9218 | 26.9035 40 | 22.5492 | 24.8648 45 | 21.657 | 23.4282 50 | 21.077 | 22.4158 55 | 20.7001 | 21.7024 60 | 20.455 | 21.1996 65 | 20.2958 | 20.8454 70 | 20.1923 | 20.5957 75 | 20.125 | 20.4198 80 | 20.0812 | 20.2958 85 | 20.0528 | 20.2085 90 | 20.0343 | 20.1469 95 | 20.0223 | 20.1035 100 | 20.0145 | 20.073 STEP: 10 Time Euler Analytic ------------------------- 0 | 100.0 | 100.0 10 | 44.0 | 59.7268 20 | 27.2 | 39.7278 30 | 22.16 | 29.7965 40 | 20.648 | 24.8648 50 | 20.1944 | 22.4158 60 | 20.0583 | 21.1996 70 | 20.0175 | 20.5957 80 | 20.0052 | 20.2958 90 | 20.0016 | 20.1469 100 | 20.0005 | 20.073
version 1
/* REXX ***************************************************************
* 24.05.2013 Walter Pachl translated from PL/I
**********************************************************************/
Numeric Digits 100
T0=100
Tr=20
k=0.07
h=2
x=t0
Call head
do t=0 to 100 by 2
Select
When t<=4 | t>=96 Then
call o x
When t=8 Then
Say '...'
Otherwise
Nop
End
x=x+h*f(x)
end
h=5
y=t0
Call head
do t=0 to 100 by 5
call o y
y=y+h*f(y)
end
h=10
z=t0
Call head
do t=0 to 100 by 10
call o z
z=z+h*f(z)
end
Exit
f: procedure Expose k Tr
Parse Arg t
return -k*(T-Tr)
head:
Say 'h='h
Say ' t By formula By Euler'
Return
o:
Parse Arg v
Say right(t,3) format(Tr+(T0-Tr)/exp(k*t),5,10) format(v,5,10)
Return
exp: Procedure
Parse Arg x,prec
If prec<9 Then prec=9
Numeric Digits (2*prec)
Numeric Fuzz 3
o=1
u=1
r=1
Do i=1 By 1
ra=r
o=o*x
u=u*i
r=r+(o/u)
If r=ra Then Leave
End
Numeric Digits (prec)
r=r+0
Return r
Output:
h=2 t By formula By Euler 0 100.0000000000 100.0000000000 2 89.5486587628 88.8000000000 4 80.4626994233 79.1680000000 ... 96 20.0965230572 20.0574137147 98 20.0839131147 20.0493757946 100 20.0729505571 20.0424631834 h=5 t By formula By Euler 0 100.0000000000 100.0000000000 5 76.3750471216 72.0000000000 10 59.7268242534 53.8000000000 15 47.9950199099 41.9700000000 20 39.7277571000 34.2805000000 25 33.9019154664 29.2823250000 30 29.7965142633 26.0335112500 35 26.9034869314 23.9217823125 40 24.8648050015 22.5491585031 45 23.4281701466 21.6569530270 50 22.4157906708 21.0770194676 55 21.7023789162 20.7000626539 60 21.1996461464 20.4550407250 65 20.8453763508 20.2957764713 70 20.5957266443 20.1922547063 75 20.4198014729 20.1249655591 80 20.2958290978 20.0812276134 85 20.2084672415 20.0527979487 90 20.1469043822 20.0343186667 95 20.1035217684 20.0223071333 100 20.0729505571 20.0144996367 h=10 t By formula By Euler 0 100.0000000000 100.0000000000 10 59.7268242534 44.0000000000 20 39.7277571000 27.2000000000 30 29.7965142633 22.1600000000 40 24.8648050015 20.6480000000 50 22.4157906708 20.1944000000 60 21.1996461464 20.0583200000 70 20.5957266443 20.0174960000 80 20.2958290978 20.0052488000 90 20.1469043822 20.0015746400 100 20.0729505571 20.0004723920
Version 2
Include How to use
Include Source code
Below program generates a series of results with increasing step size 2,5,10 and thus 50,20,10 steps.
-- 12 Jun 2026
include Setting
say "EULER METHOD FOR ODE y'=-0.07(y-20)"
say version
say
call Task 0,100,2,100
call Task 0,100,5,100
call Task 0,100,10,100
call Timer
exit
Task:
-- Euler method
procedure
arg x0,xn,dx,y0
say 'From' x0 'to' xn 'step' dx
say ' x Calc y True y Abs error Rel error'
x=x0; y=y0
do until x > xn
if x//10 = 0 then
call Show x,y
y+=dx*Dy(y); x+=dx
end
say
return
Show:
-- Display a line
procedure
arg xx,yy
yt=True(xx)
say Format(xx,3,0) Format(yy,3,7) Format(yt,3,7),
Left(Std(Abs(yy-yt)),9) Left(Std(Abs(yy/yt-1)),9)
return
Dy:
-- Differential equation
procedure
arg yy
return -0.07*(yy-20)
True:
-- Real solution
procedure
arg xx
return 20+(100-20)*Exp(-0.07*xx)
include Math
- Output:
EULER METHOD FOR ODE y'=-0.07(y-20) REXX-Regina_3.9.7(MT) 5.00 18 Mar 2025 From 0 to 100 step 2 x Calc y True y Abs error Rel error 0 100.0000000 100.0000000 0 0 10 57.6341614 59.7268243 2.0926629 0.0350372 20 37.7041263 39.7277571 2.0236308 0.0509374 30 28.3284994 29.7965143 1.4680149 0.0492680 40 23.9179512 24.8648050 0.9468538 0.0380800 50 21.8431101 22.4157907 0.5726806 0.0255480 60 20.8670488 21.1996462 0.3325974 0.0156888 70 20.4078832 20.5957266 0.1878434 0.0091205 80 20.1918793 20.2958291 0.1039498 0.0051217 90 20.0902653 20.1469044 0.0566391 0.0028113 100 20.0424634 20.0729506 0.0304872 0.0015188 From 0 to 100 step 5 x Calc y True y Abs error Rel error 0 100.0000000 100.0000000 0 0 10 53.8000000 59.7268243 5.9268243 0.0992322 20 34.2805000 39.7277571 5.4472571 0.1371146 30 26.0335113 29.7965143 3.763003 0.1262900 40 22.5491585 24.8648050 2.3156465 0.0931294 50 21.0770195 22.4157907 1.3387712 0.0597244 60 20.4550408 21.1996462 0.7446054 0.0351234 70 20.1922547 20.5957266 0.4034719 0.0195900 80 20.0812276 20.2958291 0.2146015 0.0105736 90 20.0343186 20.1469044 0.1125858 0.0055882 100 20.0144996 20.0729506 0.058451 0.0029119 From 0 to 100 step 10 x Calc y True y Abs error Rel error 0 100.0000000 100.0000000 0 0 10 44.0000000 59.7268243 15.726824 0.2633125 20 27.2000000 39.7277571 12.527757 0.3153401 30 22.1600000 29.7965143 7.6365143 0.2562888 40 20.6480000 24.8648050 4.216805 0.1695893 50 20.1944000 22.4157907 2.2213907 0.0990993 60 20.0583200 21.1996462 1.1413262 0.0538370 70 20.0174960 20.5957266 0.5782306 0.0280752 80 20.0052488 20.2958291 0.2905803 0.0143172 90 20.0015746 20.1469044 0.1453298 0.0072135 100 20.0004724 20.0729506 0.0724782 0.0036107 0.012 seconds
decimals(3)
see euler("return -0.07*(y-20)", 100, 0, 100, 2) + nl
see euler("return -0.07*(y-20)", 100, 0, 100, 5) + nl
see euler("return -0.07*(y-20)", 100, 0, 100, 10) + nl
func euler df, y, a, b, s
t = a
while t <= b
see "" + t + " " + y + nl
y += s * eval(df)
t += s
end
return yOutput:
0 100 2 88.800 4 79.168 6 70.884 8 63.761 10 57.634
This is a typical task for which RPL was designed. (tn,yn) are handled as complex numbers. This makes the iterations easier to calculate and facilitates the eventual display of the resulting curve, as RPL uses this data format to designate the pixels.
≪ → t temp '-K*(temp-TR)' ≫ ‘F’ STO
≪ → h fn
≪ 1 10 START
DUP h OVER C→R fn EVAL h * R→C +
NEXT 10 →LIST
≫ ≫ ‘EULER’ STO
0.07 'K' STO 20 'TR' STO
(0,100) 2 'F' EULER
- Output:
1: { (2,88.8) (4,79.168) (6,70.88448) (8,63.7606528) (10,57.634161408) (12,52.3653788109) (14,47.8342257774) (16,43.9374341685) (18,40.5861933849) (20,37.704126311) }
Analytical solution
≪ { } 2 20 FOR t t 'TR+(100-TR)*EXP(-K*t)' EVAL R→C + NEXT ≫
- Output:
1: { (2,89.5486588319) (3,84.8467396776) (4,80.4626993165) (5,76.3750471775) (6,72.5637455852) (7,69.0101115348) (8,65.6967251079) (9,62.6073440806) (10,59.7268243033) (11,57.0410454649) (12,54.5368418743) (13,52.2019379227) (14,50.0248879081) (15,47.9950199289) (16,46.1023835698) (17,44.3377011253) (18,42.69232212) (19,41.158180904) (20,39.7277571153) }
Accuracy to the degree for the first seconds of cooling needs h to be set at 0.2 s or below
def euler(y, a, b, h)
a.step(b,h) do |t|
puts "%7.3f %7.3f" % [t,y]
y += h * yield(t,y)
end
end
[10, 5, 2].each do |step|
puts "Step = #{step}"
euler(100,0,100,step) {|time, temp| -0.07 * (temp - 20) }
puts
end
- Output:
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 100.000 20.000 Step = 5 0.000 100.000 5.000 72.000 10.000 53.800 15.000 41.970 20.000 34.280 25.000 29.282 30.000 26.034 35.000 23.922 40.000 22.549 45.000 21.657 50.000 21.077 55.000 20.700 60.000 20.455 65.000 20.296 70.000 20.192 75.000 20.125 80.000 20.081 85.000 20.053 90.000 20.034 95.000 20.022 100.000 20.014 Step = 2 0.000 100.000 2.000 88.800 4.000 79.168 6.000 70.884 8.000 63.761 10.000 57.634 12.000 52.365 14.000 47.834 16.000 43.937 18.000 40.586 20.000 37.704 22.000 35.226 24.000 33.094 26.000 31.261 28.000 29.684 30.000 28.328 32.000 27.163 34.000 26.160 36.000 25.297 38.000 24.556 40.000 23.918 42.000 23.369 44.000 22.898 46.000 22.492 48.000 22.143 50.000 21.843 52.000 21.585 54.000 21.363 56.000 21.172 58.000 21.008 60.000 20.867 62.000 20.746 64.000 20.641 66.000 20.551 68.000 20.474 70.000 20.408 72.000 20.351 74.000 20.302 76.000 20.259 78.000 20.223 80.000 20.192 82.000 20.165 84.000 20.142 86.000 20.122 88.000 20.105 90.000 20.090 92.000 20.078 94.000 20.067 96.000 20.057 98.000 20.049 100.000 20.042
fn header() {
print!(" Time: ");
for t in (0..100).step_by(10) {
print!(" {:7}", t);
}
println!();
}
fn analytic() {
print!("Analytic: ");
for t in (0..=100).step_by(10) {
print!(" {:7.3}", 20.0 + 80.0 * (-0.07 * f64::from(t)).exp());
}
println!();
}
fn euler<F: Fn(f64) -> f64>(f: F, mut y: f64, step: usize, end: usize) {
print!(" Step {:2}: ", step);
for t in (0..=end).step_by(step) {
if t % 10 == 0 {
print!(" {:7.3}", y);
}
y += step as f64 * f(y);
}
println!();
}
fn main() {
header();
analytic();
for &i in &[2, 5, 10] {
euler(|temp| -0.07 * (temp - 20.0), 100.0, i, 100);
}
}
- Output:
Time: 0 10 20 30 40 50 60 70 80 90 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
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
See Ol. That implementation is valid Scheme.
In many Scheme implementations (Chez Scheme, Gauche Scheme, Gambit Scheme, CHICKEN Scheme if run with "-R r7rs", etc.), the program will run without modification. For some other implementations, the call to inexact must be either removed or changed to exact->inexact. (The call is necessary for ol, which otherwise would have written fractions that gnuplot does not understand.)
import <Utilities/Conversion.sl>;
import <Utilities/Sequence.sl>;
T0 := 100.0;
TR := 20.0;
k := 0.07;
main(args(2)) :=
let
results[i] := euler(newtonCooling, T0, 100, stringToInt(args[i]), 0, "delta_t = " ++ args[i]);
in
delimit(results, '\n');
newtonCooling(t) := -k * (t - TR);
euler: (float -> float) * float * int * int * int * char(1) -> char(1);
euler(f, y, n, h, x, output(1)) :=
let
newOutput := output ++ "\n\t" ++ intToString(x) ++ "\t" ++ floatToString(y, 3);
newY := y + h * f(y);
newX := x + h;
in
output when x > n
else
euler(f, newY, n, h, newX, newOutput);Based on C# version [1] but using tail recursion instead of looping.
- Output:
For step size 10:
main.exe 10
"delta_t = 10
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
100 20.000"
func euler_method(t0, t1, k, step_size) {
var results = [[0, t0]]
for s in (step_size..100 -> by(step_size)) {
t0 -= ((t0 - t1) * k * step_size)
results << [s, t0]
}
return results;
}
func analytical(t0, t1, k, time) {
(t0 - t1) * exp(-time * k) + t1
}
var (T0, T1, k) = (100, 20, .07)
var r2 = euler_method(T0, T1, k, 2).grep { _[0] %% 10 }
var r5 = euler_method(T0, T1, k, 5).grep { _[0] %% 10 }
var r10 = euler_method(T0, T1, k, 10).grep { _[0] %% 10 }
say "Time\t 2 err(%) 5 err(%) 10 err(%) Analytic"
say "-"*76
r2.range.each { |i|
var an = analytical(T0, T1, k, r2[i][0])
printf("%4d\t#{'%9.3f' * 7}\n",
r2[i][0],
r2[i][1], ( r2[i][1] / an) * 100 - 100,
r5[i][1], ( r5[i][1] / an) * 100 - 100,
r10[i][1], (r10[i][1] / an) * 100 - 100,
an)
}
- 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.281 -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
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
This program outputs commands for Gnuplot, which will produce a PNG. Run the program with a command such as poly --script name_you_gave_the_file.sml | gnuplot.
(* Approximate y(t) in dy/dt=f(t,y), y(a)=y0, t going from a to b with
positive step size h. Produce a list of point pairs as output. *)
fun eulerMethod (f, y0, a, b, h) =
let
fun loop (t, y, pointPairs) =
let
val pointPairs = (t, y) :: pointPairs
in
if b <= t then
rev pointPairs
else
loop (t + h, y + (h * f (t, y)), pointPairs)
end
in
loop (a, y0, nil)
end
(* How to step temperature according to Newton's law of cooling. *)
fun f (t, temp) = ~0.07 * (temp - 20.0)
val data2 = eulerMethod (f, 100.0, 0.0, 100.0, 2.0)
and data5 = eulerMethod (f, 100.0, 0.0, 100.0, 5.0)
and data10 = eulerMethod (f, 100.0, 0.0, 100.0, 10.0)
fun printPointPairs pointPairs =
app (fn (t, y) => (print (Real.toString t);
print " ";
print (Real.toString y);
print "\n"))
pointPairs
;
print ("set encoding utf8\n");
print ("set term png size 1000,750 font 'Brioso Pro,16'\n");
print ("set output 'newton-cooling-SML.png'\n");
print ("set grid\n");
print ("set title 'Newton\\U+2019s Law of Cooling'\n");
print ("set xlabel 'Elapsed time (seconds)'\n");
print ("set ylabel 'Temperature (Celsius)'\n");
print ("set xrange [0:100]\n");
print ("set yrange [15:100]\n");
print ("y(x) = 20.0 + (80.0 * exp (-0.07 * x))\n");
print ("plot y(x) with lines title 'Analytic solution', \\\n");
print (" '-' with linespoints title 'Euler method, step size 2s', \\\n");
print (" '-' with linespoints title 'Step size 5s', \\\n");
print (" '-' with linespoints title 'Step size 10s'\n");
printPointPairs data2;
print ("e\n");
printPointPairs data5;
print ("e\n");
printPointPairs data10;
print ("e\n");
- Output:
import Foundation
let numberFormat = " %7.3f"
let k = 0.07
let initialTemp = 100.0
let finalTemp = 20.0
let startTime = 0
let endTime = 100
func ivpEuler(function: (Double, Double) -> Double, initialValue: Double, step: Int) {
print(String(format: " Step %2d: ", step), terminator: "")
var y = initialValue
for t in stride(from: startTime, through: endTime, by: step) {
if t % 10 == 0 {
print(String(format: numberFormat, y), terminator: "")
}
y += Double(step) * function(Double(t), y)
}
print()
}
func analytic() {
print(" Time: ", terminator: "")
for t in stride(from: startTime, through: endTime, by: 10) {
print(String(format: " %7d", t), terminator: "")
}
print("\nAnalytic: ", terminator: "")
for t in stride(from: startTime, through: endTime, by: 10) {
let temp = finalTemp + (initialTemp - finalTemp) * exp(-k * Double(t))
print(String(format: numberFormat, temp), terminator: "")
}
print()
}
func cooling(t: Double, temp: Double) -> Double {
return -k * (temp - finalTemp)
}
analytic()
ivpEuler(function: cooling, initialValue: initialTemp, step: 2)
ivpEuler(function: cooling, initialValue: initialTemp, step: 5)
ivpEuler(function: cooling, initialValue: initialTemp, step: 10)
- Output:
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
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
Solution:
"Euler Solution"
T ← 100 # initial starting temp
TR ← 20 # room temp
TMINUSTR ← - TR T
h ← 10 # step size
k ← 0.07 # coefficent
TEND ← 100 # end time
n ← ÷ h 100 # steps
# inital starting point
T
.
# .. clone the top of stack and take if for next step
# repeat the steps n times with ⍥
Solution ← [⍥(.. - × h × k - TR)]+ n 1
⇌ ⊂ Solution T
# analytical solution
"Analytical Solution"
# apply function to LIST
List ← × k × h ⇡n
# Analytical solution applied
+ TR × TMINUSTR ⁿ ¯List eExample:
"Euler Solution" [100 43.99999999999999 27.199999999999996 22.159999999999997 20.648 20.194399999999998 20.05832 20.017496 20.0052488 20.00157464 20.000472392 20.0001417176 20.0001417176 20.0001417176] "Analytical Solution" [100 59.72682430331276 39.727757115328515 29.796514260238553 24.864805010017434 22.41579067378548 21.199646145638216 20.59572664567395 20.295829097318634 20.14690438216231]
Private Sub ivp_euler(f As String, y As Double, step As Integer, end_t As Integer)
Dim t As Integer
Debug.Print " Step "; step; ": ",
Do While t <= end_t
If t Mod 10 = 0 Then Debug.Print Format(y, "0.000"),
y = y + step * Application.Run(f, y)
t = t + step
Loop
Debug.Print
End Sub
Sub analytic()
Debug.Print " Time: ",
For t = 0 To 100 Step 10
Debug.Print " "; t,
Next t
Debug.Print
Debug.Print "Analytic: ",
For t = 0 To 100 Step 10
Debug.Print Format(20 + 80 * Exp(-0.07 * t), "0.000"),
Next t
Debug.Print
End Sub
Private Function cooling(temp As Double) As Double
cooling = -0.07 * (temp - 20)
End Function
Public Sub euler_method()
Dim r_cooling As String
r_cooling = "cooling"
analytic
ivp_euler r_cooling, 100, 2, 100
ivp_euler r_cooling, 100, 5, 100
ivp_euler r_cooling, 100, 10, 100
End Sub
- Output:
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,281 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
import math
// Fdy is a type for fntion f used in Euler's method.
type Fdy = fn(f64, f64) f64
// euler_step 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.
fn euler_step(f Fdy, x f64, y f64, h f64) f64 {
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.
// new_cooling_rate returns a fntion that computes cooling rate
// for a given cooling rate constant k.
fn new_cooling_rate(k f64) fn(f64) f64 {
return fn[k](delta_temp f64) f64 {
return -k * delta_temp
}
}
// new_temp_func returns a fntion that computes the analytical solution
// of cooling rate integrated over time.
fn new_temp_func(k f64, ambient_temp f64, initial_temp f64) fn(f64) f64 {
return fn[ambient_temp,initial_temp,k](time f64) f64 {
return ambient_temp + (initial_temp-ambient_temp)*math.exp(-k*time)
}
}
// new_cooling_rate_dy returns a fntion of the kind needed for Euler's method.
// That is, a fntion representing dy(x, y(x)).
//
// Parameters to new_cooling_rate_dy are cooling constant k and ambient
// temperature.
fn new_cooling_rate_dy(k f64, ambient_temp f64) Fdy {
// note that result is dependent only on the object temperature.
// there are no additional dependencies on time, so the x parameter
// provided by euler_step is unused.
return fn[k,ambient_temp](_ f64, object_temp f64) f64 {
return new_cooling_rate(k)(object_temp - ambient_temp)
}
}
fn main() {
k := .07
temp_room := 20.0
temp_object := 100.0
fcr := new_cooling_rate_dy(k, temp_room)
analytic := new_temp_func(k, temp_room, temp_object)
for delta_time in [2.0, 5, 10] {
println("Step size = ${delta_time:.1f}")
println(" Time Euler's Analytic")
mut temp := temp_object
for time := 0.0; time <= 100; time += delta_time {
println("${time:5.1f} ${temp:7.3f} ${analytic(time):7.3f}")
temp = euler_step(fcr, time, temp, delta_time)
}
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
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
import "./fmt" for Fmt
import "./iterate" for Stepped
var euler = Fn.new { |f, y, step, end|
Fmt.write(" Step $2d: ", step)
for (t in Stepped.new(0..end, step)) {
if (t%10 == 0) Fmt.write(" $7.3f", y)
y = y + step * f.call(y)
}
System.print()
}
var analytic = Fn.new {
System.write(" Time: ")
for (t in Stepped.new(0..100, 10)) Fmt.write(" $7d", t)
System.write("\nAnalytic: ")
for (t in Stepped.new(0..100, 10)) {
Fmt.write(" $7.3f", 20 + 80 * (-0.07*t).exp)
}
System.print()
}
var cooling = Fn.new { |temp| -0.07 * (temp - 20) }
analytic.call()
for (i in [2, 5, 10]) euler.call(cooling, 100, i, 100)
- Output:
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
!ys-0
defn main(*):
exact =: analytical(100.0)
steps =:: [2.0, 5.0, 10.0]
each h steps:
approx =: euler(cooling 0.0 100.0 100.0 h)
say: "h=$(h:I): $(format '%.4f' approx) (exact: $(format '%.4f' exact))"
defn analytical(t):
20.0 + (80.0 * Math/exp(-0.07 * t))
defn euler(f t0 y0 t-end h):
loop t t0, y y0:
if t >= t-end:
then: y
else:
t2 =: t + h
y2 =: y + (h * f(t y))
recur: t2 y2
defn cooling(t T):
-0.07 * (T - 20.0)
- Output:
$ ys euler-method.ys h=2: 20.0425 (exact: 20.0730) h=5: 20.0145 (exact: 20.0730) h=10: 20.0005 (exact: 20.0730)
import "std/math.zc"
alias deriv = fn*(f64, f64) -> f64;
fn euler(f: deriv, y: f64, step: int, end: int) {
let t = 0;
print " Step {step:2d}: ";
do {
if !(t % 10) { print " {y:7.3f}"; }
y += step * f(t, y);
t += step;
} while t <= end;
println "";
}
fn analytic() {
print " Time: ";
for let t = 0; t <= 100; t += 10 { print " {(f64)t:7g}"; }
print "\nAnalytic: ";
for let t = 0; t <= 100; t += 10 {
let v = 20.0 + 80.0 * Math::exp(-0.07 * (f64)t);
print " {v:7.3f}";
}
println "";
}
fn cooling(_: f64, temp: f64) -> f64 {
return -0.07 * (temp - 20.0);
}
fn main() {
analytic();
euler(cooling, 100, 2, 100);
euler(cooling, 100, 5, 100);
euler(cooling, 100, 10, 100);
}
- Output:
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
const FMT=" %7.3f";
fcn ivp_euler(f,y,step,end_t){
print(" Step %2d: ".fmt(step));
foreach t in ([0..end_t,step]){
if (t % 10 == 0) print(FMT.fmt(y));
y += f(t,y) * step;
}
println();
}
fcn analytic{
print(" Time: ");
foreach t in ([0..100,10]){ print(" %7g".fmt(t)) }
print("\nAnalytic: ");
foreach t in ([0..100,10]){ print(FMT.fmt(20.0 + 80.0 * (-0.07 * t).exp())) }
println();
}
fcn cooling(_,temp){ return(-0.07 * (temp - 20)) }
analytic();
ivp_euler(cooling, 100.0, 2, 100);
ivp_euler(cooling, 100.0, 5, 100);
ivp_euler(cooling, 100.0, 10, 100);- Output:
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
10 LET d$="-0.07*(y-20)": LET y=100: LET a=0: LET b=100: LET s=10
20 LET t=a
30 IF t<=b THEN PRINT t;TAB 10;y: LET y=y+s*VAL d$: LET t=t+s: GO TO 30Categories:
- Programming Tasks
- Mathematical operations
- 11l
- Ada
- ALGOL 68
- ALGOL W
- Arturo
- ATS
- BASIC
- Applesoft BASIC
- BASIC256
- BBC BASIC
- Chipmunk Basic
- Craft Basic
- FreeBASIC
- Gambas
- GW-BASIC
- MSX Basic
- QBasic
- Run BASIC
- True BASIC
- XBasic
- Yabasic
- C
- C sharp
- C++
- Clay
- Clojure
- COBOL
- Common Lisp
- D
- Dart
- Delphi
- DuckDB
- EasyLang
- Elixir
- Erlang
- Euler Math Toolbox
- F Sharp
- Factor
- Forth
- Fortran
- Futhark
- FutureBasic
- Go
- Groovy
- Haskell
- Icon
- Unicon
- J
- K
- Java
- JavaScript
- Jq
- Julia
- Kotlin
- Lambdatalk
- Lua
- Maple
- Mathematica
- Wolfram Language
- MATLAB
- Maxima
- МК-61/52
- Nim
- Objeck
- ObjectIcon
- OCaml
- Odin
- Oforth
- Ol
- Pascal
- PascalABC.NET
- Perl
- Phix
- Phix/pGUI
- Phix/online
- PicoLisp
- PL/I
- Pluto
- Pluto-fmt
- PowerShell
- PureBasic
- Python
- Q
- R
- Racket
- Raku
- Rebol
- REXX
- Ring
- RPL
- Ruby
- Rust
- Scala
- Scheme
- SequenceL
- Sidef
- Smalltalk
- Standard ML
- Swift
- Tcl
- Uiua
- VBA
- V (Vlang)
- XPL0
- Wren
- Wren-fmt
- Wren-iterate
- YAMLScript
- Zen C
- Zkl
- ZX Spectrum Basic
- Pages with too many expensive parser function calls