Euler method: Difference between revisions
(Ada solution added) |
|||
Line 51: | Line 51: | ||
[[Image:Euler_Method_Newton_Cooling.png|center|750px]] |
[[Image:Euler_Method_Newton_Cooling.png|center|750px]] |
||
=={{header|Ada}}== |
|||
The solution is generic, usable for any floating point type. The package specification: |
|||
<lang Ada> |
|||
generic |
|||
type Number is digits <>; |
|||
package Euler is |
|||
type Waveform is array (Integer range <>) of Number; |
|||
function Solve |
|||
( F : not null access function (T, Y : Number) return Number; |
|||
Y0 : Number; |
|||
T0, T1 : Number; |
|||
N : Positive |
|||
) return Waveform; |
|||
end Euler; |
|||
</lang> |
|||
The function Solve returns the solution of the differential equation for each of N+1 points, starting from the point T0. The implementation: |
|||
<lang Ada> |
|||
package body Euler is |
|||
function Solve |
|||
( F : not null access function (T, Y : Number) return Number; |
|||
Y0 : Number; |
|||
T0, T1 : Number; |
|||
N : Positive |
|||
) return Waveform is |
|||
dT : constant Number := (T1 - T0) / Number (N); |
|||
begin |
|||
return Y : Waveform (0..N) do |
|||
Y (0) := Y0; |
|||
for I in 1..Y'Last loop |
|||
Y (I) := Y (I - 1) + dT * F (T0 + dT * Number (I - 1), Y (I - 1)); |
|||
end loop; |
|||
end return; |
|||
end Solve; |
|||
end Euler; |
|||
</lang> |
|||
The test program: |
|||
<lang Ada> |
|||
with Ada.Text_IO; use Ada.Text_IO; |
|||
with Euler; |
|||
procedure Test_Euler_Method is |
|||
package Float_Euler is new Euler (Float); |
|||
use Float_Euler; |
|||
function Newton_Cooling_Law (T, Y : Float) return Float is |
|||
begin |
|||
return -0.07 * (Y - 20.0); |
|||
end Newton_Cooling_Law; |
|||
Y : Waveform := Solve (Newton_Cooling_Law'Access, 100.0, 0.0, 100.0, 10); |
|||
begin |
|||
for I in Y'Range loop |
|||
Put_Line (Integer'Image (10 * I) & ":" & Float'Image (Y (I))); |
|||
end loop; |
|||
end Test_Euler_Method; |
|||
</lang> |
|||
Sample output: |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|ALGOL 68}}== |
=={{header|ALGOL 68}}== |
||
{{trans|D}} Note: This specimen retains the original [[#D|D]] coding style. |
{{trans|D}} Note: This specimen retains the original [[#D|D]] coding style. |
Revision as of 14:58, 19 March 2011
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:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{dy(t)}{dt} = f(t,y(t))}
with an initial value
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle y(t_0) = y_0}
To get a numeric solution, we replace the derivative on the LHS with a finite difference approximation:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{dy(t)}{dt} \approx \frac{y(t+h)-y(t)}{h}}
then solve for Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle y(t+h)} :
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle y(t+h) \approx y(t) + h \, \frac{dy(t)}{dt}}
which is the same as
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle y(t+h) \approx y(t) + h \, f(t,y(t))}
The iterative solution rule is then:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle y_{n+1} = y_n + h \, f(t_n, y_n)}
Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle h} is the step size, the most relevant parameter for accuracy of the solution. A smaller step size increases accuracy but also the computation cost, so it has always has to be hand-picked according to the problem at hand.
Example: Newton's Cooling Law
Newton's cooling law describes how an object of initial temperature Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle T(t_0) = T_0} cools down in an environment of temperature Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle T_R } :
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{dT(t)}{dt} = -k \, \Delta T}
or
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{dT(t)}{dt} = -k \, (T(t) - T_R)}
It says that the cooling rate Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{dT(t)}{dt}} of the object is proportional to the current temperature difference Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta T = (T(t) - T_R)} to the surrounding environment.
The analytical solution, which we will compare to the numerical approximation, is
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle T(t) = T_R + (T_0 - T_R) \; e^{-k t}}
Task
The task is to implement a routine of Euler's method and then to use it to solve the given example of Newton's cooling law with it for three different step sizes of 2 s, 5 s and 10 s and to compare with the analytical solution. The initial temperature Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle T_0} shall be 100 °C, the room temperature Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle T_R} 20 °C, and the cooling constant Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} 0.07. The time interval to calculate shall be from 0 s to 100 s.
An reference solution (Common Lisp) can be seen below. We see that bigger step sizes lead to reduced approximation accuracy.
Ada
The solution is generic, usable for any floating point type. The package specification: <lang Ada> generic
type Number is digits <>;
package Euler is
type Waveform is array (Integer range <>) of Number; function Solve ( F : not null access function (T, Y : Number) return Number; Y0 : Number; T0, T1 : Number; N : Positive ) return Waveform;
end Euler; </lang> The function Solve returns the solution of the differential equation for each of N+1 points, starting from the point T0. The implementation: <lang Ada> package body Euler is
function Solve ( F : not null access function (T, Y : Number) return Number; Y0 : Number; T0, T1 : Number; N : Positive ) return Waveform is dT : constant Number := (T1 - T0) / Number (N); begin return Y : Waveform (0..N) do Y (0) := Y0; for I in 1..Y'Last loop Y (I) := Y (I - 1) + dT * F (T0 + dT * Number (I - 1), Y (I - 1)); end loop; end return; end Solve;
end Euler; </lang> The test program: <lang Ada> with Ada.Text_IO; use Ada.Text_IO; with Euler;
procedure Test_Euler_Method is
package Float_Euler is new Euler (Float); use Float_Euler;
function Newton_Cooling_Law (T, Y : Float) return Float is begin return -0.07 * (Y - 20.0); end Newton_Cooling_Law; Y : Waveform := Solve (Newton_Cooling_Law'Access, 100.0, 0.0, 100.0, 10);
begin
for I in Y'Range loop Put_Line (Integer'Image (10 * I) & ":" & Float'Image (Y (I))); end loop;
end Test_Euler_Method; </lang> Sample output:
0: 1.00000E+02 10: 4.40000E+01 20: 2.72000E+01 30: 2.21600E+01 40: 2.06480E+01 50: 2.01944E+01 60: 2.00583E+01 70: 2.00175E+01 80: 2.00052E+01 90: 2.00016E+01 100: 2.00005E+01
ALGOL 68
Note: This specimen retains the original D coding style.
<lang algol68># Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and t=a..b and the step size h.
PROC euler = (PROC(REAL,REAL)REAL f, REAL y0, a, b, h)REAL: (
REAL y := y0, t := a; WHILE t < b DO printf(($g(-6,3)": "g(-7,3)l$, t, y)); y +:= h * f(t, y); t +:= h OD; printf($"done"l$); y
);
- Example: Newton's cooling law #
PROC newton cooling law = (REAL time, t)REAL: (
-0.07 * (t - 20)
);
main: (
euler(newton cooling law, 100, 0, 100, 10)
)</lang> Ouput:
0.000: 100.000 10.000: 44.000 20.000: 27.200 30.000: 22.160 40.000: 20.648 50.000: 20.194 60.000: 20.058 70.000: 20.017 80.000: 20.005 90.000: 20.002 done
C++
<lang cpp>#include <iomanip>
- include <iostream>
typedef double F(double,double);
/* Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and t=a..b and the step size h.
- /
void euler(F f, double y0, double a, double b, double h) {
double y = y0; for (double t = a; t < b; t += h) { std::cout << std::fixed << std::setprecision(3) << t << " " << y << "\n"; y += h * f(t, y); } std::cout << "done\n";
}
// Example: Newton's cooling law double newtonCoolingLaw(double, double t) {
return -0.07 * (t - 20);
}
int main() {
euler(newtonCoolingLaw, 100, 0, 100, 2); euler(newtonCoolingLaw, 100, 0, 100, 5); euler(newtonCoolingLaw, 100, 0, 100, 10);
}</lang> Last part of output:
... 0.000 100.000 10.000 44.000 20.000 27.200 30.000 22.160 40.000 20.648 50.000 20.194 60.000 20.058 70.000 20.017 80.000 20.005 90.000 20.002 done
Common Lisp
<lang lisp>;; 't' usually means "true" in CL, but we need 't' here for time/temperature. (defconstant true 'cl:t) (shadow 't)
- Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and t=a..b and the step size h.
(defun euler (f y0 a b h)
;; Set the initial values and increments of the iteration variables. (do ((t a (incf t h)) (y y0 (incf y (* h (funcall f t y)))))
;; End the iteration when t reaches the end b of the time interval. ((>= t b) 'DONE)
;; Print t and y(t) at every step of the do loop. (format true "~6,3F ~6,3F~%" t y)))
- Example
- Newton's cooling law, f(t,T) = -0.07*(T-20)
(defun newton-cooling (time T) (* -0.07 (- T 20)))
- Generate the data for all three step sizes (2,5 and 10).
(euler #'newton-cooling 100 0 100 2) (euler #'newton-cooling 100 0 100 5) (euler #'newton-cooling 100 0 100 10)</lang>
Example output: 0.000 100.000 10.000 44.000 20.000 27.200 30.000 22.160 40.000 20.648 50.000 20.194 60.000 20.058 70.000 20.017 80.000 20.005 90.000 20.002 DONE
D
<lang d>import std.stdio, std.range;
/** Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and t=a..b and the step size h.
- /
void euler(F)(F f, double y0, double a, double b, double h) {
double y = y0; foreach (t; iota(a, b, h)) { writefln("%.3f %.3f", t, y); y += h * f(t, y); } writeln("done");
}
/// Example: Newton's cooling law auto newtonCoolingLaw(double time, double t) {
return -0.07 * (t - 20);
}
void main() {
euler(&newtonCoolingLaw, 100, 0, 100, 2); euler(&newtonCoolingLaw, 100, 0, 100, 5); euler(&newtonCoolingLaw, 100, 0, 100, 10);
}</lang> Last part of the output:
... 0.000 100.000 10.000 44.000 20.000 27.200 30.000 22.160 40.000 20.648 50.000 20.194 60.000 20.058 70.000 20.017 80.000 20.005 90.000 20.002 done
J
Solution: <lang j>NB.*euler a Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and t=a..b and step size h. euler=: adverb define
'Y a b h'=. 4{. y t=. a while. b >: tincr=. {:t + h do. Y=. Y, (+ h * u) {:Y t=. t, tincr end. t,.Y
)
ncl=: _0.07 * -&20 NB. Newton's Cooling Law </lang> Example: <lang j> ncl euler 100 0 100 2 ... NB. output redacted for brevity
ncl euler 100 0 100 5
... NB. output redacted for brevity
ncl euler 100 0 100 10 0 100 10 44 20 27.2 30 22.16 40 20.648 50 20.1944 60 20.0583 70 20.0175 80 20.0052 90 20.0016
100 20.0005</lang>
Lua
<lang lua>T0 = 100 TR = 20 k = 0.07 delta_t = { 2, 5, 10 } n = 100
NewtonCooling = function( t ) return -k * ( t - TR ) end
function Euler( f, y0, n, h )
local y = y0 for x = 0, n, h do
print( "", x, y )
y = y + h * f( y ) end
end
for i = 1, #delta_t do
print( "delta_t = ", delta_t[i] ) Euler( NewtonCooling, T0, n, delta_t[i] )
end </lang>
PicoLisp
<lang PicoLisp>(load "@lib/math.l")
(de euler (F Y A B H)
(while (> B A) (prinl (round A) " " (round Y)) (inc 'Y (*/ H (F A Y) 1.0)) (inc 'A H) ) )
(de newtonCoolingLaw (A B)
(*/ -0.07 (- B 20.) 1.0) )
(euler newtonCoolingLaw 100.0 0 100.0 2.0) (euler newtonCoolingLaw 100.0 0 100.0 5.0) (euler newtonCoolingLaw 100.0 0 100.0 10.0)</lang> Output:
... 0.000 100.000 10.000 44.000 20.000 27.200 30.000 22.160 40.000 20.648 50.000 20.194 60.000 20.058 70.000 20.018 80.000 20.005 90.000 20.002
Tcl
<lang tcl>proc euler {f y0 a b h} {
puts "computing $f over \[$a..$b\], step $h" set y [expr {double($y0)}] for {set t [expr {double($a)}]} {$t < $b} {set t [expr {$t + $h}]} {
puts [format "%.3f\t%.3f" $t $y] set y [expr {$y + $h * double([$f $t $y])}]
} puts "done"
}</lang> Demonstration with the Newton Cooling Law: <lang tcl>proc newtonCoolingLaw {time temp} {
expr {-0.07 * ($temp - 20)}
}
euler newtonCoolingLaw 100 0 100 2 euler newtonCoolingLaw 100 0 100 5 euler newtonCoolingLaw 100 0 100 10</lang> End of output:
... computing newtonCoolingLaw over [0..100], step 10 0.000 100.000 10.000 44.000 20.000 27.200 30.000 22.160 40.000 20.648 50.000 20.194 60.000 20.058 70.000 20.017 80.000 20.005 90.000 20.002 done