Solving the steady and unsteady 2D heat conduction problem
1. Objective:
Solve the 2D heat conduction equation for steady-state analysis & Transient State analysis by using point iterative techniques for both implicit and explicit schemes (Transient Analysis only). Implement the following methods
- Jacobi
- Gauss-seidel
- Successive over-relaxation
[For absolute error criteria of 1e-4]
2. Problem Definition:
- Assume that the domain is a unit square.
- Assume nx = ny [Number of points along the x-direction is equal to the number of points along the y direction]
- Boundary conditions for steady and transient case
- Left: 400K
- Top: 600K
- Right: 800K
- Bottom: 900K
3. Theory:
Explicit and implicit methods are used in numerical analysis for obtaining numerical approximations for time-dependent ordinary and partial differential equations. Heat conduction in 2D for steady and transient problem can be discretized using Taylor Series Approximation, and solved using explicit or implicit schemes.
3.a. Steady-State:
The 2D steady-state heat conduction equation is as follows:
`(∂^2 T)/(∂x^2) + (∂^2 T)/(∂y^2) =0`
Using forward differencing scheme in time and central differencing scheme in space, the above equation can be simplified as an algebraic equation to be solved using point iterartive techniques:
`(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/(∆x^2 )+ (T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/(∆y^2 )=0`
Using implicit scheme:
`T_(i,j)^(n+1)*(2/(∆x^2) + 2/(∆y^2) )= [(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))/(∆x^2) ]+ [(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))/(∆y^2) ]`
`T_(i,j)^(n+1)*(2*G1+2*G2)= G1*[T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1) ]+ G2*[T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1) ]`
On simplifying, we get the following linear equations:
`T_(i,j)^(n+1)= (G1)/(2*(G1+G2) )*[T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1) ]+( G2)/(2*(G1+G2) )*[T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1) ]`
This can be written as;
`T_(i,j)^(n+1) = β*H + ρ*V`
where,
`H= [T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1) ] ; V= [T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1) ] ; β= (G1)/(2*(G1 + G2)) ; ρ= (G2)/(2*(G1 + G2))`
3.a.1. Jacobi:
In this method values of nodal temperature in previous iteration is used to compute nodal temperature values in the current iteration.
`H= [old__T_(i+1,j)^(n+1)+old__T_(i-1,j)^(n+1)] ;V= [old__T_(i,j+1)^(n+1)+old__T_(i,j-1)^(n+1) ]`
3.a.2. Gauss Seidel:
In this method, unlike Jacobi, solution is made to converge faster by using nodal temperature values from the current iteration for (i-1,j)th and (i,j-1)th nodes when solving for (i,j)th node.
`H= [old__T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1)] ;V= [old__T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1) ]`
3.a.3. SOR:
In SOR method, convergence is achieved much faster by increasing the rate of diffusion of information within the solution domain, by multiplying with a factor alpha. Value of alpha depends on various factors.
`H= [old__T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1)] ;V= [old__T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1) ]`
`T_(i,j)^(n+1) = (1 - ALPHA)*T_(i,j)^n + ALPHA*[β*H + ρ*V]`
The Matlab code for solving the 2D steady-state heat conduction using point iterative techniques are given below. The code outputs the steady-state temperature distribution for the given problem using the methods discussed above.
clear all
close all
clc
%%%%%%%%%%%%%%%%%%%%%%% Steady State Problem %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Solving 2D Heat Conduction Problem %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assumptions:
% Unit square domain
% nx = ny
% BC: Left:400K
% Top:600K
% Right:800K
% Bottom:900K
% Absolute Error = 1e-4
%************************ Solver Parameters *************************%
Lx = 1;
Ly = 1;
nx = 81;
ny = nx;
dx = Lx/nx;
dy = Ly/ny;
tolerance = 1e-4;
%*************************** Calculations ***************************%
K1 = 1/dx^2;
K2 = 1/dy^2;
Const = 1/(2*K1 + 2*K2);
%%%%%%%%%%%%%%%%%%%%%%%%%% Implicit Schemes %%%%%%%%%%%%%%%%%%%%%%%%
for method = 1:3
error = 9e9;
%************************* Initialization ***************************%
x = linspace(0,1,nx);
y = x;
T = ones(nx,ny);
[X, Y] = meshgrid(x,y);
%*********************** Boundary Conditions ***********************%
T(1,:) = 600; % Top
T(nx,:) = 900; % Bottom
T(:,1) = 400; % Left
T(:,ny) = 800; % Right
T(1,1) = (400+600)/2; % Upper Left
T(1,nx)= (800+600)/2; % Upper Right
T(ny,1)= (400+900)/2; % Bottom Left
T(ny,nx)= (800+900)/2; % Bottom Right
T_old = T;
T_prev_dt = T;
iter = 1;
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%% a) Jacobi %%%%%%%%%%%%%%%%%%%%%%%%%%%%
if method == 1
%************************ convergence loop *************************%
while (error > tolerance)
%*********************** marching along Y *************************%
for i = 2:(ny - 1)
%*********************** marching along X *************************%
for j = 2:(nx - 1)
Horiz = T_old(i+1,j) + T_old(i-1,j);
Vert = T_old(i,j+1) + T_old(i,j-1);
T(i,j) = K1*Const*Horiz + K2*Const*Vert;
end
end
%******************** Computing maximum error *********************%
error = max(max(abs(T-T_old)));
T_old = T;
iter = iter + 1;
end
%***************** Plotting Steady State Temp. Dist. ***************%
tt = int16(toc);
s = sprintf(\'Using Jacobi, Temperature Distribution at steady state (iteration no. = %d ,time elapsed = %ds)\',iter,tt);
figure(1)
contourf(X, Y, T, \'ShowText\',\'on\')
set(gca,\'YDIR\',\'reverse\')
title(s)
grid on
axis equal
colorbar
colormap(jet)
end
%%%%%%%%%%%%%%%%%%%%%%%%% b) Gauss Seidel %%%%%%%%%%%%%%%%%%%%%%%%%
if method == 2
%************************ convergence loop *************************%
while (error > tolerance)
%*********************** marching along Y *************************%
for i = 2:(ny - 1)
%*********************** marching along X *************************%
for j = 2:(nx - 1)
Horiz = T_old(i+1,j) + T(i-1,j);
Vert = T_old(i,j+1) + T(i,j-1);
T(i,j) = K1*Const*Horiz + K2*Const*Vert;
end
end
%******************** Computing maximum error *********************%
error = max(max(abs(T-T_old)));
T_old = T;
iter = iter + 1;
end
%***************** Plotting Steady State Temp. Dist. ***************%
tt = int16(toc);
s = sprintf(\'Using Gauss Seidel, Temperature Distribution at steady state (iteration no. = %d ,time elapsed = %ds)\',iter,tt);
figure(2)
contourf(X, Y, T, \'ShowText\',\'on\')
set(gca,\'YDIR\',\'reverse\')
title(s)
grid on
axis equal
colorbar
colormap(jet)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c) SOR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha = 1.7;
if method == 3
%************************ convergence loop *************************%
while (error > tolerance)
%*********************** marching along Y *************************%
for i = 2:(ny - 1)
%*********************** marching along X *************************%
for j = 2:(nx - 1)
Horiz = T_old(i+1,j) + T(i-1,j);
Vert = T_old(i,j+1) + T(i,j-1);
T(i,j) = T_old(i,j)*(1 - alpha) + alpha*[K1*Const*Horiz + K2*Const*Vert];
end
end
%******************** Computing maximum error *********************%
error = max(max(abs(T-T_old)));
T_old = T;
iter = iter + 1;
end
%***************** Plotting Steady State Temp. Dist. ***************%
tt = int16(toc);
s = sprintf(\'Using SOR, Temperature Distribution at steady state (iteration no. = %d ,time elapsed = %ds)\',iter,tt);
figure(3)
contourf(X, Y, T, \'ShowText\',\'on\')
set(gca,\'YDIR\',\'reverse\')
title(s)
grid on
axis equal
colorbar
colormap(jet)
end
end
3.b. Transient State:
3.b.1. Explicit:
Explicit methods calculate the state of a system at a later time from the state of the system at the current time.
`(∂T)/(∂t)+ α[(∂^2 T)/(∂x^2) + (∂^2 T)/(∂y^2) ] = 0`
`(T_(i,j)^(n+1)-T_(i,j)^n)/(∆t)+ α*[(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/((∆x)^2 )+ (T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/((∆y)^2 )]=0`
`T_(i,j)^(n+1)=T_(i,j)^n- α*∆t*[(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/((∆x)^2 )+ (T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/((∆y)^2 )]^(n)`
`T_(i,j)^(n+1)=T_(i,j)^n- α*∆t*[(T_(i+1,j)^(n)-2*T_(i,j)^n+T_(i-1,j)^(n))/((∆x)^2) + (T_(i,j+1)^(n)-2*T_(i,j)^n+T_(i,j-1)^(n))/((∆y)^2) ]`
`T_(i,j)^(n+1) = T_(i,j)^n - β*H- ρ*V`
Where,
`H= [(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))/((∆x)^2) ] ;V= [(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))/((∆y)^2) ] ; β = (K1)/(1 + 2*K1 + 2*K2) ; ρ = (K2)/(1 + 2*K1 + 2*K2)`
3.b.2. Implicit:
Implicit methods find a solution by solving an equation involving both the current state of the system and the later one.
`(∂T)/(∂t)+ α[(∂^2 T)/(∂x^2) + (∂^2 T)/(∂y^2) ] = 0`
`(T_(i,j)^(n+1)-T_(i,j)^n)/(∆t)+ α*[(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/((∆x)^2 )+ (T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/((∆y)^2 )]=0`
`T_(i,j)^(n+1)=T_(i,j)^n- α*∆t*[(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/((∆x)^2 )+ (T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/((∆y)^2 )]^(n+1)`
`T_(i,j)^(n+1)*(1+ ( 2*α*∆t)/((∆x)^2) +(2*α*∆t)/((∆y)^2) )=T_(i,j)^n- α*∆t*[(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))/((∆x)^2) + (T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))/((∆y)^2) ]`
`T_(i,j)^(n+1)*(1 + 2*K1 + 2*K2)=T_(i,j)^n- K1*[(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))/((∆x)^2) ] + K2*[(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))/((∆y)^2) ]`
`T_(i,j)^(n+1) = T_(i,j)^n - (K1)/(1 + 2*K1 + 2*K2)*[(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))/((∆x)^2) ] + (K2)/(1 + 2*K1 + 2*K2)*[(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))/((∆y)^2) ]`
`T_(i,j)^(n+1) = T_(i,j)^n - β*H- ρ*V`
Where,
`H= [(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))/((∆x)^2) ] ;V= [(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))/((∆y)^2) ] ; β = (K1)/(1 + 2*K1 + 2*K2) ; ρ = (K2)/(1 + 2*K1 + 2*K2)`
3.b.2.1. Jacobi:
This formula expresses temperature at a node (i; j) in terms of its neighbouring values taken from previous iteration level. Progress of solution is monitored through absolute error for all nodes after every iteration. This an elliptic problem and solution must depend on the B.C., but this method diffuses the B.C. extremely slow.
`H= [(old__T_(i+1,j)^(n+1)+old__T_(i-1,j)^(n+1))/(∆x)^2 ] ;V= [(old__T_(i,j+1)^(n+1)+old__T_(i,j-1)^(n+1))/(∆y)^2 ]`
3.b.2.2. Gauss-Seidel:
In this method, latest values of variables are considered as soon as they are available; which causes faster rate of convergence towards solution.
`H= [(old__T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))/(∆x)^2 ] ;V= [(old__T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))/(∆y)^2 ]`
3.b.2.3. SOR:
In the SOR method, convergence is achieved much faster by increasing the rate of diffusion of information within the solution domain, by multiplying with a factor-alpha. Value of alpha depends on various factors.
`T_(i,j)^(n+1) = T_(i,j)^n - ALPHA*[β*H + ρ*V]`
The Matlab code for solving the 2D transient heat conduction using point iterative techniques are given below. The code outputs the transient temperature distribution for the specified time, for the given problem using the methods discussed above.
clear all
close all
clc
%%%%%%%%%%%%%%%%%%%%%%% Transient Problem %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Solving 2D Heat Conduction Problem %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Note:
%%%%% 1. Methods 1 to 3 are Implicit, Method 4 is Explicit
%%%%% 2. To reduce computation time, explicit solver only prints
%%%%% specified number of frames (100 by default). To print
%%%%% more frames, change the value of \"frames\" variable in the
%%%%% solver parameter.
%%%%% 3. Parameters are common for both Implicit and Explicit Solvers
%%%%% unless stated in comments.
%%%%%
% Assumptions:
% Unit square domain
% nx = ny
% BC: Left:400K
% Top:600K
% Right:800K
% Bottom:900K
% Absolute Error = 1e-4
%************************ Solver Parameters *************************%
Lx = 1;
Ly = 1;
nx = 41;
ny = nx;
%%%%% Set the no. of frames to be printed
frames = 100; %%%%%%%% For Explicit Solver
Simulation_Time = 2.5;
dt = 1e-2; %%%%%%%% For Implicit Solver
d_t = 1.4895e-4; %%%%%%%% For Explicit Solver
nt = Simulation_Time/dt; %%%%%%%% For Implicit Solver
n_t = Simulation_Time/d_t; %%%%%%%% For Explicit Solver
dx = Lx/nx;
dy = Ly/ny;
tolerance = 1e-4;
%*************************** Calculations ***************************%
%%%%%%%% For Implicit Solver
alpha = 1;
K1 = dt*alpha/dx^2;
K2 = dt*alpha/dy^2;
Const = 1/(1 + 2*K1 + 2*K2);
%%%%%%%% For Explicit Solver
G1 = d_t*alpha/dx^2;
G2 = d_t*alpha/dy^2;
marker = n_t/frames;
flag = marker;
for method = 1:4
error = 9e9;
%%%%%%%%%%%%%%%%%%%%%%%%%% Implicit Schemes %%%%%%%%%%%%%%%%%%%%%%%%
%************************* Initialization ***************************%
x = linspace(0,1,nx);
y = x;
T = ones(nx,ny);
[X, Y] = meshgrid(x,y);
%*********************** Boundary Conditions ***********************%
T(1,:) = 600; % Top
T(nx,:) = 900; % Bottom
T(:,1) = 400; % Left
T(:,ny) = 800; % Right
T(1,1) = (400+600)/2; % Upper Left
T(1,nx)= (800+600)/2; % Upper Right
T(ny,1)= (400+900)/2; % Bottom Left
T(ny,nx)= (800+900)/2; % Bottom Right
T_old = T;
T_prev_dt = T;
iter = 1;
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%% a) Jacobi %%%%%%%%%%%%%%%%%%%%%%%%%%%%
if method == 1
%**************************** Time loop ****************************%
for t = 1:nt
%************************ convergence loop *************************%
while (error > tolerance)
%*********************** marching along Y *************************%
for i = 2:(ny - 1)
%*********************** marching along X *************************%
for j = 2:(nx - 1)
Horiz = T_old(i+1,j) + T_old(i-1,j);
Vert = T_old(i,j+1) + T_old(i,j-1);
T(i,j) = Const*T_prev_dt(i,j) + K1*Const*Horiz + K2*Const*Vert;
end
end
%******************** Computing maximum error *********************%
error = max(max(abs(T-T_old)));
T_old = T;
iter = iter + 1;
end
error = 9e9;
T_prev_dt = T;
%****************** Plotting Transient Temp. Dist. ******************%
tt = int16(toc);
s = sprintf(\'Using Jacobi, Temp. Dist. at time %ds, (iteration no. = %d, time elapsed = %ds)\',dt*t,iter,tt);
figure(1)
contourf(X, Y, T, \'ShowText\',\'on\')
set(gca,\'YDIR\',\'reverse\')
title(s)
grid on
axis equal
colorbar
colormap(jet)
pause(0.01)
end
end
%%%%%%%%%%%%%%%%%%%%%%%%% b) Gauss Seidel %%%%%%%%%%%%%%%%%%%%%%%%%
if method == 2
%**************************** Time loop ****************************%
for t = 1:nt
%************************ convergence loop *************************%
while (error > tolerance)
%*********************** marching along Y *************************%
for i = 2:(ny - 1)
%*********************** marching along X *************************%
for j = 2:(nx - 1)
Horiz = T_old(i+1,j) + T(i-1,j);
Vert = T_old(i,j+1) + T(i,j-1);
T(i,j) = Const*T_prev_dt(i,j) + K1*Const*Horiz + K2*Const*Vert;
end
end
%******************** Computing maximum error *********************%
error = max(max(abs(T-T_old)));
T_old = T;
iter = iter + 1;
end
error = 9e9;
T_prev_dt = T;
%****************** Plotting Transient Temp. Dist. ******************%
tt = int16(toc);
s = sprintf(\'Using Gauss Seidel, Temp. Dist. at time %ds, (iteration no. = %d, time elapsed = %ds)\',dt*t,iter,tt);
figure(2)
contourf(X, Y, T, \'ShowText\',\'on\')
set(gca,\'YDIR\',\'reverse\')
title(s)
grid on
axis equal
colorbar
colormap(jet)
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c) SOR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha_SOR = 1.7;
if method == 3
%**************************** Time loop ****************************%
for t = 1:nt
%************************ convergence loop *************************%
while (error > tolerance)
%*********************** marching along Y *************************%
for i = 2:(ny - 1)
%*********************** marching along X *************************%
for j = 2:(nx - 1)
Horiz = T_old(i+1,j) + T(i-1,j);
Vert = T_old(i,j+1) + T(i,j-1);
T(i,j) = (Const*T_prev_dt(i,j) + K1*Const*Horiz + K2*Const*Vert)*alpha_SOR + T_old(i,j)*(1-alpha_SOR);
end
end
%******************** Computing maximum error *********************%
error = max(max(abs(T-T_old)));
T_old = T;
iter = iter + 1;
end
error = 9e9;
T_prev_dt = T;
%****************** Plotting Transient Temp. Dist. ******************%
tt = int16(toc);
s = sprintf(\'Using SOR, Temp. Dist. at time %ds, (iteration no. = %d, time elapsed = %ds)\',dt*t,iter,tt);
figure(3)
contourf(X, Y, T, \'ShowText\',\'on\')
set(gca,\'YDIR\',\'reverse\')
title(s)
grid on
axis equal
colorbar
colormap(jet)
end
end
if method == 4
%%%%%%%%%%%%%%%%%%%%%%%%%% Explicit Scheme %%%%%%%%%%%%%%%%%%%%%%%%%
%**************************** Time loop ****************************%
for t = 1:n_t
%*********************** marching along Y *************************%
for i = 2:(ny - 1)
%*********************** marching along X *************************%
for j = 2:(nx - 1)
Horiz = T_prev_dt(i+1,j) + T_prev_dt(i-1,j) - 2*T_prev_dt(i,j);
Vert = T_prev_dt(i,j+1) + T_prev_dt(i,j-1) - 2*T_prev_dt(i,j);
T(i,j) = T_prev_dt(i,j) + G1*Horiz + G2*Vert;
end
end
T_prev_dt = T;
%****************** Plotting Transient Temp. Dist. ******************%
if (t >= flag)||(int32(t) == int32(n_t))
tt = int16(toc);
s = sprintf(\' Temperature Distribution at time %ds (time elapsed = %ds)\',d_t*t,tt);
figure(4)
contourf(X, Y, T, \'ShowText\',\'on\')
set(gca,\'YDIR\',\'reverse\')
title(s)
grid on
axis equal
colorbar
colormap(jet)
flag = flag + marker;
end
end
end
end
4. Result and Discussions:
4.a. Steady-state
4.a.1. Jacobi
- Number of iterations: 11721
- Least preferable method due to a large number of iteration required for achieving steady-state
- Diffusion of B.C. is slow
4.a.2. Gauss seidel
- Number of iterations: 6317
- Significantly faster than Jacobi
- Faster diffusion of B.C.
4.a.3. SOR
- Number of iterations: 1309
- The most preferable method for steady-state simulation
- Diffusion of B.C. is accelerated by using a factor \'alpha\'
- Diffusion rate can be varied by varying \'alpha\'
4.b. Transient:
- Simulation time: 2.5s
4.b.1. Explicit:
- The explicit scheme requires a \'dt\' value of at least 1.4895e-4 to obtain a stable solution
- This method is not prefered due to its small time-step value which increases computation time
- Plotting every time step in the 0-2.5s range takes more than 80 minutes, hence to reduce the time required for simulation only 100 equally spaced time steps in the 0-2.5s range is being plotted. This reduces the elapsed time to be much lower than the implicit method. This is because the implicit method plots every time step in the 0-2.5s range, hence time elapsed is not a suitable parameter to evaluate the performance of implicit and explicit methods.
4.b.2. Implicit:
- Faster than the explicit method, and can obtain reasonably accurate values at higher time steps
4.b.2.1. Jacobi:
- Number of iterations: 17146
- Slowest implicit method
4.b.2.2. Gauss-Seidel:
- Number of iterations: 9884
- This method is significantly faster than Jacobi
4.b.2.3. SOR:
- Number of iterations: 2478
- The most preferable method for transient simulation
5. Conclusions:
- It can be concluded that for the steady-state and transient problem, SOR is the best method for solving this problem as it requires significantly fewer iterations compared to both Jacobi and Gauss-Seidel methods (Note: The time elapsed is not an appropriate measure of performance as inferred from the result in the explicit method. The number of iterations is a suitable measure of solver performance in this case)
- Even though the explicit method doesn\'t have a convergence loop, any gain in simulation time is offset by the smaller size of time step required to obtain a reasonably accurate result
- To reduce the computation time of the explicit method, 1.4895e-4s is the largest time step value that can be used to a stable solution for a simulation time of 5s