MATLAB Modelling of Heat Equation

April 29, 2022

By Shubham

Table of Contents

1

2

The Heat Equation: Mathematical Preliminaries

Obtaining the heat equation:

Initial and boundary conditions:

Solving the heat equation:

3

Modelling the heat equation in MATLAB:

4

Analyzing the result:

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.

Be a Brand Representative

Join Training Modules

Participate in Online Quiz

Get Research Assistance

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.

Get instant access to the code, model, or application of the video or article you found helpful! Simply purchase the specific title, if available, and receive the download link right away! #MATLABHelper #CodeMadeEasy

Ready to take your MATLAB skills to the next level? Look no further! At MATLAB Helper, we've got you covered. From free community support to expert help and training, we've got all the resources you need to become a pro in no time. If you have any questions or queries, don't hesitate to reach out to us. Simply post a comment below or send us an email at [email protected].

And don't forget to connect with us on LinkedIn, Facebook, and Subscribe to our YouTube Channel! We're always sharing helpful tips and updates, so you can stay up-to-date on everything related to MATLAB. Plus, if you spot any bugs or errors on our website, just let us know and we'll make sure to fix it ASAP.

Ready to get started? Book your expert help with Research Assistance plan today and get personalized assistance tailored to your needs. Or, if you're looking for more comprehensive training, join one of our training modules and get hands-on experience with the latest techniques and technologies. The choice is yours – start learning and growing with MATLAB Helper today!

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.

Need Urgent Help?

Our experts assist in all MATLAB & Simulink fields with communication options from live sessions to offline work.

Book Expert Help

testimonials

Philippa E. / PhD Fellow

I STRONGLY recommend MATLAB Helper to EVERYONE interested in doing a successful project & research work! MATLAB Helper has completely surpassed my expectations. Just book their service and forget all your worries.

Yogesh Mangal / Graduate Trainee

MATLAB Helper provide training and internship in MATLAB. It covers many topics of MATLAB. I have received my training from MATLAB Helper with the best experience. It also provide many webinar which is helpful to learning in MATLAB.

Research Assistance

Training Programs

Online Courses

USEFUL resources