.st0{fill:#FFFFFF;}

MATLAB Modelling of Heat Equation 

 April 29, 2022

By  Shubham

Introduction:

According to Wikipedia, a partial differential equation (PDE) is an equation that imposes relations between the various partial derivatives of a multivariable function. PDE forms the foundational basis for modern scientific understanding of sound, heat, diffusion, thermodynamics, and quantum mechanics. This blog will deal with applying partial differential equations in the form of the heat equation and its modelling in MATLAB.

The Heat Equation: Mathematical Preliminaries

Obtaining the heat equation:

Let’s start with the derivation of the heat equation. Consider a heat-conducting homogeneous rod, extending from x=0 to x=L along the x-axis having a uniform cross-section A.
metallic rod

Heat conduction in a thin bar

Let the material density be \rho . Assume the rod is laterally insulated so that heat flows only in the direction of the x-axis. Let u(x,t) denote the temperature of the cross-section at point x at any instant of time t and sigma be the specific heat of the rod.

The amount of heat in the rod between the cross-section at x and x+\Delta x is,

Q(t)=\int_{x}^{x+\Delta x}\sigma \rho Au(s,t)

Rate of heat flowing into the segment x is,

-KA\frac{\partial u(x,t)}{\partial x}

Rate of heat flowing out of the segment x+\Delta x is,

-KA\frac{\partial u(x+\Delta x,t)}{\partial x}

 Here K denotes the thermal conductivity of the rod. You might be wondering about the negative sign in front of the two equations given above. The negative sign indicates that the heat flows in the direction of colder regions. The amount of heat content between the section x and x+\Delta x must be equal to the difference between the amount of heat that flows in at cross-section x and the amount heat that flows out the at cross-section x+\Delta x. Hence,

\frac{\partial Q}{\partial x}=\int_{x}^{x+\Delta x}\sigma \rho A\frac{\partial u(s,t)}{\partial x}ds

=KA\left [ \frac{\partial u(x+\Delta x,t)}{\partial x} -\frac{\partial u(x,t)}{\partial x}\right ]

To solve the integration, assume the function to be a continuous function of s and apply the mean value theorem.

\int_{x}^{x+\Delta x}\frac{\partial u(s,t)}{\partial t}ds=\frac{\partial u(\xi ,t)}{\partial x}\Delta x, x<\xi <x+\Delta x

\sigma\rho \Delta x\frac{\partial u(\xi ,t)}{\partial t}=K\left [ \frac{\partial u(x+\Delta x,t)}{\partial x} -\frac{\partial u(x,t)}{\partial x}\right ]

\frac{\partial u(x,t)}{\partial t}=a^2\frac{\partial^2 u(x,t)}{\partial x^2}

 This equation is called one-dimensional heat equation. Here, a^2 =K/(\sigma \rho ). This constant a^2 is called the diffusivity within the solid.

Initial and boundary conditions:

An initial boundary value problem for the heat equation consists of the PDE along with three other conditions specified at x=0, x=L and t=0. More formally it can be written as:

\frac{\partial u}{\partial x}=a^2\frac{\partial^2 u}{\partial x^2},0<x<L,0<t

and it describes the temperature distribution in the rod at any time 0<t subject to conditions

u(x,0)=f(x),0<x<L(i)

And

u(0,t)=u(L,t)=0,0<t(ii)

Equation (i) describes the initial conditions while equation (ii) describes the boundary conditions.

Solving the heat equation:

let’s reconsider our heat equation. We have to find the solution to the homogeneous heat equation

\frac{\partial u}{\partial t}=a^2 \frac{\partial^2 u}{\partial x^2},0<x<L,0<t

Having the initial condition u(x,0)=f(x),0<x<L and boundary conditions u(0,t)=u(L,t)=0,0<t.

We will use the variable separation method to solve the equation. The solution that we seek will be in the form u(x,t)=X(x)T(t)

Since, \frac{\partial u}{\partial t}=X(x)T'(t) and \frac{\partial^2 u}{\partial x^2}=X''T(t) . Therefore, the above equation becomes T'(t)X(x)=a^2 X''(x)T(t). This equation can be rewritten as \frac{T'}{a^2T}=\frac{X''}{X}=-\lambda . Here, -\lambda is a constant. This equation gives us two ordinary differential equations: X''+\lambda X=0 and T'+a^2\lambda T=0. The boundary condition changes to u(0,t)=X(0)T(t) and u(L,t)=X(L)T(t)=0 for 0<t.

Depending on the value of \lambda , there are three cases: \lambda=-m^2 , \lambda=0 and \lambda=k^2 . For all the three cases the general solution for the first differential equation will be as follows:

  • For \lambda=-m^2 , X(x)=Acosh(mx)+Bsinh(mx)
  • For \lambda=0 , X(x)=C+D(x)
  • For \lambda=k^2 , X(x)=Ecos(kx)+Fsin(kx)

For second differential equation, we will use \lambda _{n}=n^2\pi ^2/L^2 in the equation T'+a^2\lambda T=0. The general solution is T_{n}(t)=G_{n}exp\left ( -\frac{a^2n^2\pi ^2}{L^2}t \right ).

The particular solution of the considered heat equation is u(x,t)=\sum_{n=1}^{\infty }Csin\left ( \frac{n \pi x}{L}exp\left ( -\frac{a^2n^2\pi^2}{L^2}t \right ) \right ). The general solutions hence equal a linear sum of the particular solution. Here the coefficient C can be obtained by setting t=0 for the entire equation. Hence,

f(x)=\sum_{n=1}^{\infty }Csin\left ( \frac{n\pi x}{L} \right ),0<x<L

And,

C=\frac{2}{L}\int_{0}^{L}f(x)sin\left ( \frac{n \pi x}{L} \right )dx, n=1,2,3,...

Now that the heat equation is ready let’s model it in MATLAB.

The initialization conditions and boundary values are based on the heat equation we used and derived earlier. Hence, if the boundary conditions are different, they should be changed in the codes. With this info in mind, let’s dive into the MATLAB code.

Modelling the heat equation in MATLAB:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Shubham Kumar, MATLAB Helper
% Topic: Heat Equation
% Website: https://matlabhelper.com
% Date: 03-04-2022
% Version: R2021b
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;
clear;
close all;

% Initialisation
M = 15;
dx = 0.125;
dt = 0.05;
X = 0:dx:pi;
T = 0:dt:2;
u = zeros(length(T),length(X));

% Boundary Conditions
Xx = repmat(X,[length(T) 1]);
Tt = repmat(T',[1 length(X)]);

% Solving the equations
for m = 1:M
temp = 2*m-1;
coeff = 8 / (pi * temp * temp * temp);
u = u + coeff * sin(temp*Xx) .* exp(-temp * temp * Tt);
end

% Plotting the heat equation
surf(Xx,Tt,u)
xlabel('distance','Fontsize',20);
ylabel('time','Fontsize',20);
zlabel('u(x,t)','Fontsize',20);

Seeing the code above, one might wonder what the repmat() function is and why it is used. As per the MathWorks website, the repmet() function returns copies of an array. Hence, we are making copies of the u(x,t) function for the different values of X and T. To solve the equations, we will introduce a for loop which will go up to the value of M and calculates the coefficients and hence the function u(x,t).

Finally, let’s plot the equations and check our result.

On running the above code, MATLAB will generate the following graph:

heat equation graph

Graph of heat equation

Analyzing the result:

Since the result is a 3D plot, it can be rotated to a different point of view and analyzed. The unrotated plot tells us that temperature u(x,t) within a thin bar is zero at the ends. This result satisfies the boundary conditions. Now, let’s rotate the graph.

heat equation result rotated horizontally

Graph of heat equation rotated horizontally

The graph shows that the heat flows out from the center of the thin metallic bar towards the colder region and the ends of the metallic bar where the temperature is removed.

heat equation result rotated sidewise

Graph of heat equation rotated sidewise

The parabolic shape is reduced as we move away from the hotter region towards the colder area with the increase in time. This again shows that the heat is being removed from the ends of the metallic bar.

Did you find some helpful content from our video or article and now looking for its code, model, or application? You can purchase the specific Title, if available, and instantly get the download link.

Thank you for reading this blog. Do share this blog if you found it helpful. If you have any queries, post them in the comments or contact us by emailing your questions to [email protected]. Follow us on LinkedIn Facebook, and Subscribe to our YouTube Channel. If you find any bug or error on this or any other page on our website, please inform us & we will correct it.

If you are looking for free help, you can post your comment below & wait for any community member to respond, which is not guaranteed. You can book Expert Help, a paid service, and get assistance in your requirement. If your timeline allows, we recommend you book the Research Assistance plan. If you want to get trained in MATLAB or Simulink, you may join one of our training modules. 

If you are ready for the paid service, share your requirement with necessary attachments & inform us about any Service preference along with the timeline. Once evaluated, we will revert to you with more details and the next suggested step.

Education is our future. MATLAB is our feature. Happy MATLABing!

About the author 

Shubham

I have done my bachelor degree in electronics from University of Delhi. I like to tinker with the hardware to understand how it actually works.

{"email":"Email address invalid","url":"Website address invalid","required":"Required field missing"}

Connect with MATLAB Helper ®

Follow: YouTube Channel, LinkedIn Company, Facebook Page, Instagram Page

Join Community of MATLAB Enthusiasts: Facebook Group, Telegram, LinkedIn Group

Use Website Chat or WhatsApp at +91-8104622179

Watch our Latest Video

Meet Industrial Experts & Ask your Questions

Industrial Interaction MATLAB Helper
>