Courses by Software
Courses by Semester
Courses by Domain
Tool-focused Courses
Machine learning
POPULAR COURSES
Post Graduate Program in Hybrid Electric Vehicle Design and AnalysisPost Graduate Program in Computational Fluid DynamicsPost Graduate Program in CADPost Graduate Program in CAEPost Graduate Program in Manufacturing DesignPost Graduate Program in Computational Design and Pre-processingPost Graduate Program in Complete Passenger Car Design & Product DevelopmentSuccess Stories
MATLAB Scripting of steady and unsteady 2D heat conduction equation using Jacobi, Gauss-Seidel & SOR Method
Objective Write a MATLAB Script to solver 2D Heat Diffusion Equation for Steady-state & Transient State using Jacobi, Gauss-seidel & Successive over-relaxation iterative method for Steady-state & Implicit and Explicit schemes for Transient state. Input 2D, Plate with negligible thickness Length of Plate…
Rajat Walia
updated on 01 Feb 2021
Objective
Write a MATLAB Script to solver 2D Heat Diffusion Equation for Steady-state & Transient State using Jacobi, Gauss-seidel & Successive over-relaxation iterative method for Steady-state & Implicit and Explicit schemes for Transient state.
Input
2D, Plate with negligible thickness
Length of Plate = 1 meter, Width of Plate = 1 meter.
Convegence Criteria is 1e-4.
Boundary Condition:-
Top Boundary = 600 K, Bottom Boundary = 900 K, Left Boundary = 400 K, Right Boundary = 800 K
Uniform grid spacing along x & y direction with 10 nodes.
Thermal Diffusivity - 1.4 m2/s
Over-relaxation factor - 1.1
Timestep - 1000 (for transient simulation)
Discretization of 2D Diffusion Equation (Steady Form)
d2Tdx2+d2Tdy2=0d2Tdx2+d2Tdy2=0
Discretization scheme for diffusion term: Central Difference Scheme
Ti+1,j-2Ti,j+Ti-1,j(δx)2+Ti,j+1-2Ti,j+Ti,j-1(δy)2
Iterarive method such as Jacobi, Gauss-seidel & Successive over-relaxation can be used to solve the above equation to obtain the steady state temperature field.
Jacobi Method - Calculates the temperature of nodes at the current iteration using temperature values of the previous iteration.
Ti,j=(1k)(Toldi+1,j+Toldi-1,j(δx)2+Toldi,j+1+Toldi,j-1(δy)2)
T - Temperature of the current iteration
Told - Temperature of the previous iteration
K=2(δx2+δy2δx2δy2)
Gauss-Seidel Method - Calculates the temperature of nodes at current iteration using temperature values of the latest updated iteration during the iterative procedures.
Ti,j=(1k)(Toldi+1,j+Ti-1,j(δx)2+Toldi,j+1+Ti,j-1(δy)2)
T - Temperature of the current iteration
Told - Temperature of the previous iteration
K=2(δx2+δy2δx2δy2)
Successive over-relaxation Method - Calculates the temperature of nodes at current iteration using temperature values of the previous iteration + Rate of change of temperature with respect to iteration multiplied with an over-relaxation factor.
Ti,j=Toldi,j+α(dTdi)
where (dTdi)=TGS-Told(n+1)-n
α - over-relaxation factor
n+1 - current iteration
n - old iteration
Ti,j - Temperature at current iteration
Toldi,j - Temperature at old iteration
TGS - Temperature calculated at current iteration using Gauss-Seidel Method
Finally substituting dTdi in main formula and doing some re-arrangement, The final form we will get as follow
Ti,j=Toldi,j(1-α)+αTGS
Discretization of 2D Diffusion Equation (Transient Form)
dTdt=α(d2Tdx2+d2Tdy2)
α - Thermal Diffusivity
Explicit Method - Temperature of nodes at current time step is calculated based on previous time step nodes temperature.
(Ti,j)n+1-(Ti,j)nδt=α((Ti+1,j)n-2(Ti,j)n+(Ti-1,j)n(δx)2+(Ti,j+1)n-2(Ti,j)n+(Ti,j-1)n(δy)2)
Right hand side, all values are known so it can be directly solved for next timestep using time marching approach.
Using Jacobi Method for space marching, the final discretized equation is as follow:-
(Ti,j)n+1=(Ti,j)n+αδt(δx)2((Ti+1,j)n-2(Ti,j)n+(Ti-1,j)n)+αδt(δy)2((Ti,j+1)n-2(Ti,j)n+(Ti,j-1)n)
(Ti,j)n - Temperature at previous time step
(Ti,j)n+1 - Temperature at current time step
Stability criteria:-
αδt(δx)2+αδt(δy)2≤0.5 OR (CFLx+CFLy≤0.5)
Implicit Method - Temperature of nodes at the current time step is calculated based on the current time step nodes temperature itself. This leads to a system of linear algebirc equations.
(Ti,j)n+1-(Ti,j)nδt=α((Ti+1,j)n+1-2(Ti,j)n+1+(Ti-1,j)n+1(δx)2+(Ti,j+1)n+1-2(Ti,j)n+1+(Ti,j-1)n+1(δy)2)
Right hand side all values are unknown, Hence Iterative technique such as Jacobi, Gauss-seidel & Successive over-relaxation has to be used followed by time marching.
Using different iterative methods(Jacobi, Gauss-Seidel, SOR), Solution will end up with a different number of iterations to reach specified solution time(0.5 seconds in this case) due to the iterative process nature.
Jacobi Method
(Ti,j)n+1=(Ti,j)n+αδt(δx)2((Ti+1,j)n+1-2(Ti,j)n+1+(Ti-1,j)n+1)+αδt(δy)2((Ti,j+1)n+1-2(Ti,j)n+1+(Ti,j-1)n+1)
αδt(δx)2-CFLx
αδt(δy)2-CFLy
(Ti,j)n+1=(Ti,j)n+CFLx((Ti+1,j)n+1-2(Ti,j)n+1+(Ti-1,j)n+1)+CFLy((Ti,j+1)n+1-2(Ti,j)n+1+(Ti,j-1)n+1)
(Ti,j)n+1(1+2CFLx+2CFLy)=(Ti,j)n+CFLx((Ti+1,j)n+1+(Ti-1,j)n+1)+CFLy((Ti,j+1)n+1+(Ti,j-1)n+1)
(Ti,j)n+1=1(1+2CFLx+2CFLy)((Ti,j)n+CFLx((Ti+1,j)n+1+(Ti-1,j)n+1)+CFLy((Ti,j+1)n+1+(Ti,j-1)n+1))
1(1+2CFLx+2CFLy)-term1
CFLx⋅term1-term2
CFLy⋅term1-term3
We can write final equation for jacobi method as:-
(Ti,j)n+1=term1(Tprevi,j)n+term2((Toldi+1,j)n+1+(Toldi-1,j)n+1)+term3((Toldi,j+1)n+1+(Toldi,j-1)n+1))
(Ti,j)n+1 - Temperature at current time step
(Tprevi,j)n - Temperature at previous time step(This will remain constant with in an iterative loop and will be updated for next time loop)
(Toldi,j)n+1 - Temperature at previous iteration(This will change with in an iterative loop)
For the Gauss-Seidel method above equation can be used but with the latest updated temperature values as shown below
(Ti,j)n+1=term1(Tprevi,j)n+term2((Toldi+1,j)n+1+(Ti-1,j)n+1)+term3((Toldi,j+1)n+1+(Ti,j-1)n+1))
For the SOR method, We will use equation as follow:-
Ti,j=Toldi,j(1-α)+αTGS
TGS-Temperature calculated at current iteration using Gauss-Seidel Method
α - over-relaxation factor
MATLAB Script
I have written separate code for different type of simulations i.e. Steady simulation with Jacobi Method, Steady simulation with Gauss-Seidal Method, Steady simulation with SOR Method, Transient simulation using Explicit Method, Transient simulation using Implicit Method.
I have written a main code where I have defined all the inputs and I have scripted user defined input to select the type of simulation steady or unsteady, type of solver Jacobi, Gauss-Seidal, SOR, Explicit or Implicit.
Main Code
clc
clear all
close all
L = 1; %length of plate
W = 1; %width of plate
nx = 10; %mesh points along x direction
ny = 10; %mesh points along y direction
nt = 500; %number of time step for transient simulation
dt = 0.001; %time step size
%generating the grid
x = linspace(0,L,nx);
y = linspace(0,W,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
error = 5e9; %initial error
convergence_criteria = 1e-4; %convergence criteria
T = 300*ones(nx,ny); %initializing the temeperature field
%defining boundary condition
T_top = 600;
T_bottom = 900;
T_right = 800;
T_left = 400;
T(:,1) = T_left;
T(:,end) = T_right;
T(1,:) = T_top;
T(end,:) = T_bottom;
%corner nodes temperature
T(1,1) = (T_top + T_left)/2;
T(end,1) = (T_bottom + T_left)/2;
T(1,end) = (T_top + T_right)/2;
T(end,end) = (T_bottom + T_right)/2;
%Selectthe solver
type = sprintf("Steady Simulation - 1 nTransient Simulation - 2n");
disp(type)
simulation_type = input("Select the simulation type: 1 or 2: ");
if simulation_type == 1 %steady
solver = sprintf("nJacobi Solver - 1 nGauss-Seidel Solver - 2nSOR Solver - 3n");
disp(solver)
solver_type = input("Select the solver: 1 or 2 or 3: ");
if solver_type == 1 %steady jacobi
jacobi_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
elseif solver_type == 2 %steady gauss-seidel
gauss_seidel_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
elseif solver_type == 3 %steady SOR
SOR_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
end
elseif simulation_type == 2 %transient
solver = sprintf("nExplicit Solver - 1 nImplicit Solver - 2n");
disp(solver)
solver_type = input("Select the solver: 1 or 2: ");
if solver_type == 1 %explicit method
jacobi_unsteady_explicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria)
elseif solver_type == 2 %implicit method
sub_solver = sprintf("nJacobi Solver - 1 nGauss-Seidel Solver - 2nSOR Solver - 3n");
disp(sub_solver)
sub_solver_type = input("Select the sub solver: 1 or 2 or 3: ");
if sub_solver_type == 1 %implicit jacobi method
jacobi_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
elseif sub_solver_type == 2 %implicit gauss-seidel method
gauss_seidel_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
elseif sub_solver_type == 3 %implicit SOR method
SOR_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
end
end
endCode for Steady-State Jacobi Method
function out = jacobi_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
T_old = T; %temperature of previous iteration
jacobi_iteration = 1; %initial iteration
k = (2*(dx^2+dy^2))/((dx^2)*(dy^2));
%solving
while (max(error)>convergence_criteria)
for i = 2:nx-1 %march along x direction
for j = 2:ny-1 %march along y direction
T(i,j) = (1/k)*((T_old(i+1,j)+T_old(i-1,j))/(dx^2) + (T_old(i,j+1)+T_old(i,j-1))/(dy^2));
end
end
error = (max(abs(T_old - T))); %calculating the error
T_old = T; %updating temperature field
jacobi_iteration = jacobi_iteration + 1; %updating the iteration no.
end
%Ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
title_text = sprintf('2D Steady Heat Diffusion, Jacobi Solver, Iterations = %d',jacobi_iteration);
title(title_text);
iter = sprintf("Number of Iterations to converge: %d",jacobi_iteration);
disp(iter)
set(gca,'ydir','reverse')
endCode for Steady-State Gauss-Seidel Method
function out = gauss_seidel_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
T_old = T; %temperature of previous iteration
gauss_seidel_iteration = 1; %initial iteration
k = (2*(dx^2+dy^2))/((dx^2)*(dy^2));
%solving
while (max(error)>convergence_criteria)
for i = 2:nx-1 %march along x direction
for j = 2:ny-1 %march along y direction
T(i,j) = (1/k)*((T_old(i+1,j)+T(i-1,j))/(dx^2) + (T_old(i,j+1)+T(i,j-1))/(dy^2));
end
end
error = (max(abs(T_old - T))); %calculating the error
T_old = T; %updating temperature field
gauss_seidel_iteration = gauss_seidel_iteration + 1; %updating the iteration no.
end
%Ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
title_text = sprintf('2D Steady Heat Diffusion, Gauss-Seidel Solver, Iterations = %d',gauss_seidel_iteration);
title(title_text);
iter = sprintf("Number of Iterations to converge: %d",gauss_seidel_iteration);
disp(iter)
set(gca,'ydir','reverse')
endCode for Steady-State SOR Method
function out = SOR_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
T_old = T; %temperature of previous iteration
SOR_iteration = 1; %initial iteration
k = (2*(dx^2+dy^2))/((dx^2)*(dy^2));
alpha = 1.6; %Over-relaxation factor
%solving
while (max(error)>convergence_criteria)
for i = 2:nx-1 %march along x direction
for j = 2:ny-1 %march along y direction
T(i,j) = T(i,j)*(1-alpha) + (alpha*((1/k)*((T_old(i+1,j)+T(i-1,j))/(dx^2) + (T_old(i,j+1)+T(i,j-1))/(dy^2))));
end
end
error = (max(abs(T_old - T))); %calculating the error
T_old = T; %updating temperature field
SOR_iteration = SOR_iteration + 1; %updating iteration no.
end
%Ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
title_text = sprintf('2D Steady Heat Diffusion, SOR Solver, Iterations = %d',SOR_iteration);
title(title_text);
iter = sprintf("Number of Iterations to converge: %d",SOR_iteration);
disp(iter)
set(gca,'ydir','reverse')
endCode for Transient Explicit Method
function out = jacobi_unsteady_explicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria)
T_old = T; %temperature of previous time step
alpha = 1.4; %thermal diffusivity
CFL_x = ((2*alpha*dt)/dx^2); %CFL along x
CFL_y = ((2*alpha*dt)/dy^2); %CFL along y
solution_time = dt; %tracking the solution time
for k = 1:nt %time marching
for i = 2:nx-1 %space marching along x
for j = 2:ny-1 %space marching along y
T(i,j) = T_old(i,j) + CFL_x*(T_old(i+1,j)-2*T_old(i,j)+T_old(i-1,j)) + CFL_y*(T_old(i,j+1)-2*T_old(i,j)+T_old(i,j-1));
end
end
error = max(max(abs(T_old - T))); %calculating the error
T_old = T; %updating the temperature
%ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
solution_time = dt + solution_time;
title_text = sprintf('2D Explicit Unsteady Heat Diffusion, Jacobi, Solution Time: %.3f s',solution_time);
title(title_text);
set(gca,'ydir','reverse')
if (error < convergence_criteria);
disp("Simulation has reached steady state!")
break
end
endCode for Transient Implicit Jacobi Method
function out = jacobi_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
T_old = T; %previous iteration temperature
T_prev = T; %previous time step temperature
alpha = 1.4; %thermal diffusivity
CFL_x = ((2*alpha*dt)/dx^2); %CFL along x
CFL_y = ((2*alpha*dt)/dy^2); %CFL along y
term_1 = (1 + 2*CFL_x + 2*CFL_y)^(-1);
term_2 = term_1*CFL_x;
term_3 = term_1*CFL_y;
solution_time = dt; %tracking the solution time
jacobi_iteration = 1;
for k = 1:nt %time marching
while (error > convergence_criteria) %Jacobi Iterative Method
for i = 2:nx-1 %space marching along x
for j = 2:ny-1 %space marching along y
H = (T_old(i+1,j) + T_old(i-1,j));
V = (T_old(i,j+1) + T_old(i,j-1));
T(i,j) = ((T_prev(i,j)*term_1) + (term_2*H) + (term_3*V));
end
end
error = max(max(abs(T_old - T))); %calculating the error
T_old = T; %updating the temperature for next iteration
jacobi_iteration = 1 + jacobi_iteration;
end
error = 1; %re-assigning the error
T_prev = T; %updating the temperature for the next time step
%ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
solution_time = dt + solution_time;
title_text = sprintf('2D Implicit Unsteady Heat Diffusion, Jacobi, Solution Time: %.3f s',solution_time);
title(title_text);
set(gca,'ydir','reverse')
end
endCode for Transient Implicit Gauss-Seidel Method
function out = gauss_seidel_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
T_old = T; %previous iteration temperature
T_prev = T; %previous time step temperature
alpha = 1.4; %thermal diffusivity
CFL_x = ((2*alpha*dt)/dx^2); %CFL along x
CFL_y = ((2*alpha*dt)/dy^2); %CFL along y
term_1 = (1 + 2*CFL_x + 2*CFL_y)^(-1);
term_2 = term_1*CFL_x;
term_3 = term_1*CFL_y;
solution_time = dt; %tracking the solution time
gauss_seidel_iteration = 1;
for k = 1:nt %time marching
while (error > convergence_criteria) %Gauss-Seidel Iterative Method
for i = 2:nx-1 %space marching along x
for j = 2:ny-1 %space marching along y
H = (T_old(i+1,j) + T(i-1,j));
V = (T_old(i,j+1) + T(i,j-1));
T(i,j) = ((T_prev(i,j)*term_1) + (term_2*H) + (term_3*V));
end
end
error = max(max(abs(T_old - T))); %calculating the error
T_old = T; %updating the temperature for next iteration
gauss_seidel_iteration = 1 + gauss_seidel_iteration;
end
error = 1; %re-assigning the error
T_prev = T; %updating the temperature for the next time step
%ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
solution_time = dt + solution_time;
title_text = sprintf('2D Implicit Unsteady Heat Diffusion, Gauss-Seidel, Solution Time: %.3f s',solution_time);
title(title_text);
set(gca,'ydir','reverse')
end
endCode for Transient Implicit SOR Method
function out = SOR_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
T_old = T; %previous iteration temperature
T_prev = T; %previous time step temperature value
alpha = 1.4; %thermal diffusivity
omega = 1.1; %over-relaxation factor
CFL_x = ((2*alpha*dt)/dx^2); %CFL along x
CFL_y = ((2*alpha*dt)/dy^2); %CFL along y
term_1 = (1 + 2*CFL_x + 2*CFL_y)^(-1);
term_2 = term_1*CFL_x;
term_3 = term_1*CFL_y;
solution_time = dt; %tracking the solution time
SOR_iteration = 1;
for k = 1:nt %time marching
while (error > convergence_criteria) %SOR Iterative Method
for i = 2:nx-1 %space marching along x
for j = 2:ny-1 %space marching along y
H = (T_old(i+1,j) + T(i-1,j));
V = (T_old(i,j+1) + T(i,j-1));
T_gauss_seidel = ((T_prev(i,j)*term_1) + (term_2*H) + (term_3*V));
T(i,j) = T_old(i,j)*(1-omega) + (T_gauss_seidel*omega);
end
end
error = max(max(abs(T_old - T))); %calculating the error
T_old = T; %updating the temperature for next iteration
SOR_iteration = 1 + SOR_iteration;
end
error = 1; %re-assigning the error
T_prev = T; %updating the temperature for the next time step
%ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
solution_time = dt + solution_time;
title_text = sprintf('2D Implicit Unsteady Heat Diffusion, SOR, Solution Time: %.3f s',solution_time);
title(title_text);
set(gca,'ydir','reverse')
end
endResults for Steady Cases
| Solver Type | Number of iterations to reach convergence |
| Jacobi | 208 |
| Gauss-Seidel | 112 |
| SOR (over relexation factor = 1.1) | 38 |
SOR takes the least iterations to reach a steady-state solution whereas Jacobi takes the maximum number of iterations.
Results for Transient Cases
Explicit Scheme with Jacobi Solver
The temperature field reaches a steady-state in 0.228 seconds.
Implicit Scheme with Jacobi, Gauss-Seidel & SOR Solver
For the implicit method, the solution time is 0.5 seconds.
The temperature field has reached its steady state.
Jacobi method takes a maximum of 2318 iterations to reach 0.5 seconds of solution time.
Gauss-Seidel method takes 1790 iterations to reach 0.5 seconds of solution time.
SOR method takes a minimum of 1548 iterations to reach 0.5 seconds of solution time.
Leave a comment
Thanks for choosing to leave a comment. Please keep in mind that all the comments are moderated as per our comment policy, and your email will not be published for privacy reasons. Please leave a personal & meaningful conversation.
Other comments...
Be the first to add a comment
Read more Projects by Rajat Walia (10)
Week 9 - FVM Literature Review
Finite Volume Method Finite Volume Method is a discretization method of discretizing fluid flow governing equations. This method involves the volume integration of the governing fluid flow equation over a finite control volume. The complete domain is divided into small finite control volume cells. Governing equations…
12 Jun 2021 06:56 PM IST
Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Objective To solve Quasi 1D Converging-Diverging Nozzle numerically for both conservative & non-conservative forms of governing equations by implementing MacCormack's Method using MATLAB. Problem Description We consider a steady, isentropic flow through the converging-diverging nozzle. The inlet of the…
16 May 2021 12:00 PM IST
MATLAB Scripting of steady and unsteady 2D heat conduction equation using Jacobi, Gauss-Seidel & SOR Method
Objective Write a MATLAB Script to solver 2D Heat Diffusion Equation for Steady-state & Transient State using Jacobi, Gauss-seidel & Successive over-relaxation iterative method for Steady-state & Implicit and Explicit schemes for Transient state. Input 2D, Plate with negligible thickness Length of Plate…
01 Feb 2021 05:47 PM IST
MATLAB Scripting & Stability analysis of Linear Convection Equation
Objective Write a MATLAB Code solving Linear Convection Equation using Finite Difference Method. Study the stability of numerical scheme and analyzing the numerical diffusion error. Input 1D Domain Length - 1 m Constant Velocity(C) - 1 m/s Total Solution Time - 0.4 seconds Time Step - 0.01 seconds Number of Grid Points …
23 Jan 2021 04:36 PM IST
Related Courses
Skill-Lync offers industry relevant advanced engineering courses for engineering students by partnering with industry experts.
Our Company
4th Floor, BLOCK-B, Velachery - Tambaram Main Rd, Ram Nagar South, Madipakkam, Chennai, Tamil Nadu 600042.
Top Individual Courses
Top PG Programs
Skill-Lync Plus
Trending Blogs
© 2025 Skill-Lync Inc. All Rights Reserved.