39 lines
739 B
Plaintext
39 lines
739 B
Plaintext
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"]);
|