Euler method: Difference between revisions
(→{{header|Common Lisp}}: Added some comments.) |
|||
Line 76: | Line 76: | ||
=={{header|ALGOL 68}}== |
|||
{{trans|D}} Note: This specimen retains the original [[#D|D]] coding style. |
|||
{{works with|ALGOL 68|Revision 1 - no extensions to language used.}} |
|||
{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-1.18.0/algol68g-1.18.0-9h.tiny.el5.centos.fc11.i386.rpm/download 1.18.0-9h.tiny].}} |
|||
{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of '''format'''[ted] ''transput''.}} |
|||
<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: |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Common Lisp}}== |
=={{header|Common Lisp}}== |
||
<lang lisp> |
<lang lisp> |
Revision as of 22:02, 6 March 2011
You are encouraged to solve this task according to the task description, using any language you may know.
Method description
Euler's method numerically approximates solutions of first-order ordinary differential equations (ODEs) with a given initial value. It is a explicit method for solving initial value problems (IVPs), as described in Euler method
The ODE has to be provided in the folloving 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:
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
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 shall be 100 °C, the room temperature 20 °C, and the cooling constant 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.
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
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
--Avi 12:33, 6 March 2011 (UTC)
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