As we saw earlier, Newton's Law of Cooling states[br][br][quote]The temperature of a body changes at a rate proportional to the difference in temperature between the body and its surroundings[/quote]We saw in an earlier lesson that the above principle can be expressed as the differential equation[br][br][math]\frac{dT}{dt}=-k\left(T-T_S\right)[/math][br][br]Last time we encountered this differential equation, we used Euler's Numerical Method to get an estimate for the solution of this differential equation with [math]k=ln(5/4)[/math], an initial condition of [math]T(0)=-18[/math] Celsius, and an ambient temperature [math]T_S=+32[/math] degrees Celsius. [br][br]We will now see how to solve this differential equation algebraically as a linear differential equation. We will work with Newton's Law of Cooling in its general form, meaning we'll leave [i]k[/i] and[math]T_S[/math] unspecified, and we will [b]not [/b]have an initial condition. After we get the general solution of a Newton's Law of Cooling model, we'll then see how to work that into a real world application. [br][br]Before we can study the differential equation as a linear differential equation, we must put the differential equation in standard (linear) form. Go back and review standard linear form if you need to. [br][br][math]\frac{dT}{dt}=-k\left(T-T_S\right)[/math][br][br][math]\frac{dT}{dt}=-k\cdot T+k\cdot T_S[/math][br][br][math]\frac{dT}{dt}+k\cdot T=k\cdot T_S[/math][br][br]With the differential equation in standard (linear) form, we see that [math]p(t)=k[/math] and [math]g(t)=k\cdot T_S[/math]. In other words, both [i]p[/i] and [i]g[/i] are constant functions. As such, this will proceed more or less like the example from the earlier lesson on linear differential equation. The difference here however, is that instead of specific numbers in the equation, here we leave the constants [i]k[/i] and [math]T_S[/math] unspecified.[br][br][list=1][*][math]\int p\left(t\right)dt=\int k\cdot dt=k\cdot t[/math][br][/*][*][math]\mu\left(t\right)=e^{\int p\left(t\right)dt}=e^{k\cdot t}[/math][br][/*][*][math]T\left(t\right)=\frac{\int e^{\int p\left(t\right)dt}g\left(t\right)dt+c}{e^{\int p\left(t\right)dt}}=\frac{\int e^{k\cdot t}\cdot k\cdot T_S\cdot dt+c}{e^{k\cdot t}}[/math][br][br][br][/*][*] [math]\int e^{\int p\left(t\right)dt}g\left(t\right)dt=\int e^{k\cdot t}\cdot k\cdot T_S\cdot dt=\frac{k\cdot T_S}{k}e^{k\cdot t}=T_Se^{k\cdot t}[/math][br][/*][*][math]T\left(t\right)=\frac{\int e^{k\cdot t}\cdot k\cdot T_S\cdot dt+c}{e^{k\cdot t}}=\frac{T_Se^{k\cdot t}+c}{e^{k\cdot t}}=T_S+\frac{c}{e^{k\cdot t}}=T_S+c\cdot e^{-k\cdot t}[/math][br][/*][*]In its general form, Newton's Law of Cooling has no initial condition, and so no specific solution is required. [/*][/list][br]In practice, in a Newton's Law of Cooling situation, you are usually given the initial temperature of the body, the temperature of the body at sometime in the future, and the temperature of the surroundings ([math]T_S[/math]). You are usually [b]not [/b]provided[i] k[/i]. It turns out that these three pieces of information are just enough information to determine k and c. See the example below. [br][br][br]

[math][/math]Let's return to the example we saw earlier:[br][br]A tub of ice cream is removed from the freezer at -18 degrees Celsius and brought into a room which is (positive) 32 degrees Celsius. The temperature of the ice cream, after one minute is observed to have warmed to -8 degrees Celsius. Find the specific algebraic solution that adheres to Newton's Law of Cooling, and which models this tub of ice cream. Use the model to predict the temperature two minutes after the tub of ice cream has been moved into the 32 degree air. [br][br]SOLUTION:[br][br]First of all, let's translate the above into mathematical notation. In this setup we're looking for a function [math]T(t)[/math] which models the temperature of the ice cream, and adheres to Newton's Law of Cooling. That means that [math]T(0)=-18[/math], [math]T(1)=-8[/math] and [br][br][math]\frac{dT}{dt}=k\left(T-32\right)[/math][br][br]Unlike last time we encountered this question, we actually do not know [i]k[/i] yet. It turns out that the fact that we know the value of [i]T[/i] at two points in time ([i]t[/i]=0 and [i]t[/i]=1) is enough information to find [i]k.[/i][br][br]Now, we need not repeat the steps we did above to find the algebraic solution to the equation in standard (linear) form. Referring to step 5 from above, we know that [br][br][math]T\left(t\right)=32+ce^{-k\cdot t}[/math][br][br]solves the differential equation.[br][br]Incorporating the initial condition, [math]T(0)=-18[/math], produces the following[br][br][math]T\left(0\right)=-18=32+c\cdot e^{-k\cdot0}=32+c[/math][br][br]which can be solved for c, implying that c = - 50. [br][br]Incorporating the temperature of the ice cream after 1 minute, [math]T(1)=-8[/math], produces[br][br][math]T\left(1\right)=-8=32-50\cdot e^{-k\cdot1}=32-50\cdot e^{-k}[/math][br][br]In other words[br][br][math]-8=32-50\cdot e^{-k}[/math][br][br]We can solve this for k [br][br][math]-40=-50\cdot e^{-k}[/math][br][br][math]\frac{4}{5}=e^{-k}[/math][br][br][math]\ln\left(\frac{4}{5}\right)=-k[/math][br][br][math]k=-\ln\left(\frac{4}{5}\right)[/math][br][br]We therefore have all the information to assemble the final model of the tub of ice cream[br][br][math]T\left(t\right)=32-50e^{\ln\left(\frac{4}{5}\right)t}[/math][br][br]Plugging in [math]t=2[/math] produces the predicted temperature at minute 2, and it turns out to be 0 (check it yourself!).[br][br]Let's examine this specific algebraic solution in Geogebra, and compare it to the Euler's Numerical Method estimate of the solution we produced in a previous class. Not bad! Notice that Euler's Method is systematically over estimating the algebraic solution because the shape of the slope field points to a solution which is concave down.

Now we switch our attention back to Numerical Methods. In particular, we'll take a closer look at Euler's Method and see how to implement it programmatically so that it doesn't represent as much work as doing the calculations by hand. [br][br]For the purposes of this lesson, the most important thing to recall is the Euler's Method formulas for obtaining the coordinates of the next estimate point. [br][br][math]x_{n+1}=x_n+h[/math] [br][br][math]y_{n+1}=y_n+f\left(x_n,y_n\right)\cdot h[/math]

In this activity we will see how to implement Euler's (Numerical) Method programmatically in both Geogebra and Octave (an open source version of MATLAB). [br][br]First let's take a look at a non-linear and non-separable first order differential equation with initial condition. Feel free to use the Pen Tool to sketch an estimate of the specific solution satisfying the initial condition. Also, take a moment to reflect on why this is neither separable nor linear.[br][br][math]\frac{dy}{dx}=\sqrt{x^2+y^2};y\left(0\right)=1[/math][br][br]I've loaded this differential equation into Geogebra below. Note that I have specified the right hand side of the equation as [code]f(x,y)=sqrt(x^2+y^2)[/code], and then used [i]f[/i] when creating the slope field with [code]slopefield(f)[/code]. This is slightly easier for a number of reasons, but mainly it makes implementation of Euler's Method easier below. [br][br]

We'll implement Euler's Method in the Spreadsheet Mode of Geogebra. If you are using Geogebra Desktop, click the "View" drop down menu, and select Spreadsheet to make this mode visible. I've made it visible below.

The spreadsheet mode of Geogebra acts like a spreadsheet in Excel or Google Drive, except you can access and engage with mathematical objects in the Graphics and Algebra panes in each cell of the spreadsheet. In some sense, you can think of the Spreadsheet like a lot of input bars. [br][br]For instance, in cells [code]A1[/code] and [code]B1[/code] I've loaded the [i]x[/i] and [i]y[/i] coordinates of Initial condition with [code]=x(InitialCondition[/code]) and [code]=y(InitialCondition)[/code]

Now we'll implement Euler's Method formulas in cells [i][code]A2[/code][/i], [i][code]B2[/code][/i]. First we need a step size. We'll use [i]h[/i]=0.1 for a step size. I created this by simply typing [code]h=0.1[/code] in the input bar. [br][br]After doing this, in cell [code]A2[/code], I've calculated the [i]x[/i] coordinate of the next estimate point with[br][br][code]=A1+h[/code][br][br]And in in cell [code]B2[/code] I've calculated the [i]y[/i] coordinate of the next estimate point with [br][br][code]=B1+f(A1,B1)*h[/code][br][br]Note that because [i]f[/i] was already declared, it is easier to code the next [i]y[/i] coordinate. This is a general principle in any computer programming task: declare processes as much as possible so you don't have retype something over and over again. [br][br]Lastly I've plotted the next estimate point by entering the following in cell [code]C3[/code][br][br][code]=(A2,B2)[/code][br][br]Try out calculating the next [i]x[/i] and [i]y[/i] coordinates, and the next estimate point by putting some code in cells [code]A3[/code], [code]B3[/code] and [code]C3[/code]. Be sure to start with an equals sign![br][br]That said, don't spend too much time since we'll see how to copy the code from [code]A2[/code] and [code]B2[/code] programmatically and quickly in the next step.[br][br]

Here's another look at the same applet with the spreadsheet editor closed.

Now comes the most powerful step. Rather than typing this equation in over and over again, GeoGebra lets us apply the formula in cells [code]A2[/code], [code]B2[/code] and [code]C2[/code] for as many steps as we want by using the "fill handle" (Google it to learn more!). In most spreadsheet software (including GeoGebra and Excel), the "fill handle" is in the lower right corner of a cell, and when you drag it down, it "fills" the same formula down as far as you drag it. I've applied the fill handle to calculate 10 total estimate points. In GeoGebra it's a small blue square in the lower right of a cell, or a set of highlighted cells.[br][br]Try calculating 10 more rows by dragging the fill handle down on columns [code]A[/code], [code]B[/code], and [code]C[/code].

And again, here's a view of the same applet without the spreadsheet window open.

Not a bad fit! [br][br]Now let's see how to do this in Octave. There's nothing special about Octave; any programming language can implement Euler's Method in a "[url=https://en.wikipedia.org/wiki/For_loop]for-loop[/url]". [br][br]Click here: [url=https://octave-online.net/bucket~LxkNgcm11TXPdGU8GNCtoR]https://octave-online.net/bucket~LxkNgcm11TXPdGU8GNCtoR [/url]

You should now be able to replicate this process in Geogebra Desktop. You will need to change [code]f(x,y)[/code], the step size [i]h[/i], and also the number of steps in the spreadsheet mode, but this process should be extensible to a wide variety of use cases. That said, for your benefit, I've whipped up the following Euler's Method Calculator[br][br][url=https://www.geogebra.org/m/d6cryyfx]https://www.geogebra.org/m/d6cryyfx[/url][br][br]Looking ahead, when we encounter the Runge Kutta numerical method later this semester, we'll see that it is very similar to Euler's Method. The difference is that instead of a single formula for the next y-coordinate, we'll need a chain of several formulas to calculate the next y-coordinate. We'll discuss this later, but for now know that Runge Kutta makes use of some theory surrounding polynomial approximations to curve interpolation to get a better estimate of what the next y-coordinate should be rather than simply following the slope field in a straight line. To be continued.