Newton’s Law of Cooling with Euler’s Method Solution to IVP ODE

Topic: Newton’s Law of Cooling and Euler’s Method to IVP ODE

Subject: Physics and Numerical Methods

Tool: Scilab

By: Gani Comia, Jan. 2025

This article presents a solution as to how long it will take for an object to cool down to a specified temperature and thereby to ambient temperature. The object is being taken out of the oven at \(230^\circ C\) and put in place where the ambient temperature is \(30^\circ C\). Neglecting the heat transfer due to conduction where the object is in contact with, the time to cool it down to a particular temperature due to natural convection and radiation can be estimated using the ordinary differential equation (ODE) representing the Newton’s Law of Cooling shown in Equation (1).

$$\frac{dT}{dt}\;=\;-k \left(T(t) - T_{env} \right) \tag{1}$$

Where:

\(T(t)\) - object temperature at time \(t\),

\(T_{env}\) - environment temperature,

\(k\) - cooling constant, 

\(dT/dt\) - rate of change of temperature with time,

The analytical solution to the ODE representing Newton’s law of cooling can be found by integrating both sides of the equation as shown in Equation (2).

$$\int_{T_{0}}^{T_{t}} dT\;=\;-k \int_{0}^{t} \left(T(t)-T_{env} \right) dt \tag{2}$$

For this article, a numerical solution using the Euler method of approximating the solution of a first-order differential equation is presented with the use of Scilab. Euler methods require an ODE to be in the form of Equation (3).

$$\frac{dT}{dt}\;=\;f\left(t,T\right) \tag{3}$$

With the given initial condition,

$$\left(t_{0},\;T_{0}\right) \tag{4}$$

The approximate solution using Euler's method is shown in Equation (5).

$$\frac{\Delta T}{\Delta t}\; \approx \; \frac{T_{n+1}-T_{n}}{h}\; \approx \; f\left(t,T \right) \tag{5}$$

Or by rearranging and solving for \(T_{n+1}\),

$$T_{n+1}\;=\;T_{n}+h \;f\left(t,T\right) \tag{6}$$

The \( f(t,T)\) requires finding the cooling constant, \(k\), for a specific material. If this property seems to be unavailable or difficult to find, experimental data can be used as the basis for calculating the \(k\) using Equation (7).

$$k\;=\;- \frac{1}{t} \ln \left[ \frac{T(t) - T_{env}}{T_{0} - T_{env}} \right] \tag{7}$$

Figure (1) is the illustration of the condition to be given a solution using Newton’s law of cooling by applying Euler's method of solving ODE.

Figure 1. An object at \(230^\circ C\) exposed in an environment at \(30^\circ C\).

To provide a simulation of the given condition in Figure (1), experimental data were gathered to calculate the cooling constant, \(k\), from the given temperature and time. The cooling constant, \(k = 0.08 \;\text{per min}\), was calculated from the given parameters and Scilab script.

// Copyright (C) 2024 - Gani Comia
// Experimental Cooling Constant
clear;clc;

// (1) defining cooling constant function
function k=f(t, Tt, T0, Tenv)
    // where:
    // t - min, time consideration
    // Tt - deg C, temperature at time t
    // T0 - deg C, initial temperature
    // Tenv - deg C, environment temperature
    k = -(1/t).*log((Tt - Tenv)./(T0 - Tenv));
endfunction

// (2) experimental data
t = 30;                     // min, time between T0 and Tt
T0 = 230;                   // deg C, initial temp
Tt = 50;                    // deg C, temp at time t
Tenv = 30;                  // deg C, ambient temp

// (3) cooling constant calculation
k = f(t,Tt,T0,Tenv)
mprintf("Cooling Constant, k = %4.2f per min\n", k)

Having this cooling constant, a simulation of Newton’s law of cooling can be performed. The codes provided in the simulation should be saved on a different file from the script used for estimating the cooling constant.

The first block of the script defines the physical parameter under consideration.

clear;clc;

// (1) physical parameters
T0 = 230;                   // deg C, initial temp of the object
Tenv = 30;                  // deg C, ambient temp
k = 0.08;                   // per min, experimental cooling constant
dt = 0.05;                  // min, time step
totalTime = 80;             // min, total simulation time

The time domain in the simulation is divided into smaller step sizes for accuracy of solution. Initial time and temperature are given on the next block of script.

// (2) initial conditions and derived parameters
time = 0:dt:totalTime;
nSteps = length(time);
T = zeros(nSteps, 1);
T(1) = T0;

Equation (6), which is the numerical solution to the ODE using the Euler method, is transformed in the Scilab script. It uses the \(for-loop\) function to obtain the solution over the desired range.

// (3) numerical solution to ODE using Euler Method
for i = 1:nSteps-1
    dTdt = -k * (T(i) - Tenv);  // deg C/min, rate of temp change
    T(i+1) = T(i) + dTdt * dt;  // deg C, temp at the next time step
end

The Scilab \(plot()\) function is used to visualize the solution from the Euler method of approximation. Finding the temperature for a given time can be estimated using interpolation by the use of the \(interp1()\) function.

// (4) graph for visualization of solution
clf;
plot(time, T, "b","linewidth",2.5);
xlabel("Time, t [ min ]","fontsize",2);
ylabel("Temperature, T [ °C ]","fontsize",2);
title("Newton''s Law of Cooling Simulation","fontsize",3);
note = ["T0="+string(T0)+"°C, Tenv="+string(Tenv)+"°C, k="+string(k)]
legend(note,with_box = %F);
xgrid(3,1);

// (5) interpolating temp, Tx, for the particular time, tx
format(6)
tx = 17.5;                      // min, sample time
Tx = interp1(time, T, tx);      // deg C, interpolated temp
disp(tx, Tx)

// (6) plotting lines for tx and Tx
xptVer = [tx tx];
yptVer = [0 Tx];
plot(xptVer, yptVer, "r--","linewidth",1.5)
xptHor = [0 tx];
yptHor = [Tx Tx];
plot(xptHor, yptHor, "r--","linewidth",1.5)
plot(tx,Tx,"ro");
note1 = ["t (min) ="+string(tx)+"; T (°C) ="+string(Tx)];
xstring(tx,Tx,note1);
xstring(45,200,["https://gani-mech-toolbox.blogspot.com"])
eqn1=["$\LARGE\frac{dT}{dt}\;=\;{-k}\left(T-T_\mathbf{{env}}\right)$"];
xstring(45,166,eqn1);
a = gca();
a.children.children.mark_background = 9; 

Figure 2 is the solution over the range of time.

Figure 2. Plot of the Solution to Newton’s Law of Cooling using Euler’s Method.

This Scilab toolbox is a powerful tool for visualizing solution data and calculating specific conditions.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.