58 lines
1.7 KiB
Plaintext
58 lines
1.7 KiB
Plaintext
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
|