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