by R. Grothmann
This is just a fun application. A Pana cotta (could be any other pudding) is cooling down. The measurements were
at 0 min : 66°
at 30 min : 53.5°
at 60 min : 44.6°
The measurements are somewhat inaccurate, since the Pana cotta was stirred, windows were opened etc. However, what is the room temperature, and when will the Pana cotta have reached 35°?
We use the following model for our cooling.
whera a is the temperature of the room.
>function f(t,a,b,c) := a+exp(b*t-c)
Enter the measurements.
>t=[0,30,60]; T=[66,53.5,44.6]; >plot2d(t,T,points=1); plot2d(t,T,add=1):
Now we want to interpolate the function. We use Nelder-Mead for this.
The following is the error.
>function delta (v) := totalmax(abs(T-f(t,v[1],v[2],v[3])))
It turns out, that we need a good starting point. So we use a guess for a, the room temperature, and fit b and c according to log(M-a) = bt-c.
>a=20; M=t'|-1; b=log(T-a)'; >w=fit(M,b); b=w[1], c=w[2]
-0.0104315825592 -3.82725856953
This does not quite interpolate.
>delta([a,b,c])
0.0927776426688
So we apply Nelder-Mead. The room temperature should be 22.5°.
>w=neldermin("delta",[a,b,c]), a=w[1]; b=w[2]; c=w[3];
[22.5972, -0.0113226, -3.77052]
Plot the result.
>plot2d("f(x,a,b,c)",0,240); plot2d(t,T,points=1,add=1):
Solve for the time. After 156 min the Pana cotta should have 30°.
>solve("f(x,a,b,c)",120,y=30)
156.207181302
The true room temperature was about 24°. Let us try to fix that with another measurement.
>t=t|85; T=T|40.2; >w=neldermin("delta",[a,b,c]), a=w[1]; b=w[2]; c=w[3];
[28.2825, -0.0138396, -3.63589]
So we get another room temperature. Probably the method of measurement should have been more scientific.
>plot2d("f(x,a,b,c)",0,240); plot2d(t,T,points=1,add=1):
The predicted time is very sensitive to these errors, since it is close to the room temperature.
>solve("f(x,a,b,c)",120,y=30)
223.637373862
After 125 minutes, the predicted temperature is 34.6°.
>f(125,a,b,c)
35.0083063303
In reality, it was 35°. Let us add this to our data.
>t=t|125; T=T|35.0; >w=neldermin("delta",[a,b,c]), a=w[1]; b=w[2]; c=w[3];
[28.4844, -0.0139632, -3.63019]
>plot2d("f(x,a,b,c)",0,240); plot2d(t,T,points=1,add=1):
>solve("f(x,a,b,c)",120,y=30)
230.201564176