Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

                               Solving Steady and Unsteady State 2d Heat Conduction Equation

The two dimension heat conduction equation can be written as

`((delT)/(delt))+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`

In the above equation T represents Temperature, 'alpha' represents thermal diffusivity. x and y are spatial coordinates and t is the temporal coordinate. 

There are two approaches to solve the problem. These are listed below.

1) Steady State Equation:-

In this approach the change of temperature with time is zero, i.e the term `((delT)/(delt))' is equal to zero. The equation reduces to

`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`

In order to discretize the above equation we create a domain of unit length in x and y direction resulting in a unit square. This domain is divided 50 equally spaced grid points in the x and y direction. The Temperature boundary conditions are applied as given below

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

 

A second order central differencing scheme is applied to calculate the second order derivative terms in the steady state equation. The value of Temperature at point P is evaluated using the neighbourhood values, i.e TL and TR in x direction, Tt and Tb in the y direction.

The discretised equation for 2D steady state is given below

`(T_L-2T_p+T_R)/(Deltax^2) +(T_T-2T_p+T_B)/(Deltay^2)=0`

Solving for 'T_p'

`T_p=((Deltax^2*Deltay^2)/(2(Deltax^2+Deltay^2)))*(((T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2)))`

The term  '((2(Deltax^2+Deltay^2))/(Deltax^2*Deltay^2)))'  is labelled as k

The equation Reduces to

`T_p=(1/K)*(((T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2)))`

There are three iterative solvers which have been used to solve the above equation

 a) Jacobi Method

 In this method, given a current approximation for x, 

`x^k=(x_1^k,x_2^k,x_3^k,.......x_n^k)`

the strategy of Jacobi's Method is to use the current values of  'x_2^k','x_3^k','x_n^k' to find a new value for x_1^(k+1)

The Code written for iterative solver using Jacobi Method is given below

%Jacobi Method

clear all
close all
clc

%Input Parameter
nx=50;
ny=50;
x=linspace(0,1,nx);%Grid along X
y=linspace(0,1,ny);%Grid along y
dx=x(2)-x(1);
dy=y(2)-y(1);
T=ones(nx,ny);

%Boundary Conditions
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
Tinitial=T;
Told=T;
k=2*((dx^2+dy^2)/(dx^2*dy^2));

Tol=1e-4
error1=9e9
error2=9e9
error3=9e9
jacobi_iter=1;

tic
a=1;
while (error1>Tol)
    for i=2:nx-1
        for j=2:ny-1
            m=(Told(i+1,j)+Told(i-1,j))/(dx*dx);
            n=(Told(i,j+1)+Told(i,j-1))/(dy*dy);
            T(i,j)=(m+n)/k ;
        end
    end
    error1= max(max(abs(Told-T)));
    ERROR1(a)=error1; %Plotting error vs iterations
    Told=T;
    jacobi_iter=jacobi_iter+1;
    JI(a)=jacobi_iter;%Storing iterations in counter
    time1=toc;
    a=a+1;
    

end
figure(1);
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit1=sprintf('Number of iteration=%d\n Simulation Time=%d s',jacobi_iter,time1)
title({'Jacobi solver for steady state',tit1})
colorbar

pause(0.3)

 

-Plot for Temperature Field . For a grid spacing of 50x50, dx and dy value of 0.0204

No of iterations= 4873. Simulation time 2.189e-1seconds

b) Gauss-Seidel Method

After we compute 'x_1^(k+1)' and it is a better approximation to the true value of x_1 than x_1^k we use the new value 'x_1^(k+1)' in finding 'x_2^(k+1)' . We progressively use the updated values of the variables in the equation. The result is we get faster convergence and lesser simulation time.

Code for Gauss Seidel( TO be read in continuation of Jacobi code)

%Gauss Seidel

tic;
T=Tinitial;
Told=Tinitial;
gauss_seidel=1;
b=1;
while (error2>Tol)
    for i=2:nx-1
        for j=2:ny-1
            m=(T(i+1,j)+T(i-1,j))/(dx*dx);
            n=(T(i,j+1)+T(i,j-1))/(dy*dy);
            T(i,j)=(m+n)/k ;
        end
    end
    error2= max(max(abs(Told-T)));
      Told=T;
      ERROR2(b)= error2;
       GS(b)=gauss_seidel;
    gauss_seidel=gauss_seidel+1;
   b=b+1;
    time2=toc;

end
figure(2)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',gauss_seidel,time2)
title({'Gauss Seidel solver for steady state',tit2})
colorbar
pause(0.3)

-Plot for Temperature Field . For a grid spacing of 50x50, dx and dy value of 0.0204

No of iterations= 2610. Simulation time 1.789e-1seconds

-Successive Over Relaxation Method (SOR)

The SOR method of solving a linear system of equations Ax=B is derived by extrapolating the Gauss Seidel Method. The extrapolation takes a the form of weighted average between previous iterate and the computed Gauss seidel iterate successively for each component

`x_i^k=omega*x_i^k+(1-omega)*x_i^(k-1)`

'omega' is the relaxation factor

'omega>1', the equation is overelaxed

'omega=1', The equation is same as Gauss Seidel

'omega<1', The equation is under relaxed

 

COde for SOR (TO be read in continuation of code for Gauss Seidel)

%Successive Over Relaxation
SOR=1;

T=Tinitial
Told=Tinitial;
omega=1.5;
tic

c=1;
while (error3>Tol)
    for i=2:nx-1
        for j=2:ny-1
            m=(T(i+1,j)+T(i-1,j))/(dx*dx);
            n=(T(i,j+1)+T(i,j-1))/(dy*dy);
            
            T(i,j)=(1-omega)*Told(i,j) +omega*((m+n)/k );
        end
    end
    error3= max(max(abs(Told-T)));
    ERROR3(c)= error3;
    Told=T;
    SR(c)=SOR;
    SOR=SOR+1;
    time3=toc;
    c=c+1;

end
figure(3)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit3=sprintf('Number of iteration=%d\n Simulation Time=%d s',SOR,time3)
title({'SOR solver for steady state',tit3})
colorbar
pause(0.3)

Plot for Temperature Field . For a grid spacing of 50x50, dx and dy value of 0.0204, omega=1.5

No of iterations= 957. Simulation time 8.3477e-02seconds

 

- Plot showing variation of error with iterations for the three methods

By observing the below plot we note that SOR method converges  at much faster rate compared to Gauss seidel and Jacobi.

Unsteady Transient Equation

The unsteady Transient equation for heat conduction is given as

`(delT)/(delt) =alpha*((del^2T)/(delx^2) +(del^2T)/(dely^2) )`

Discretizing this equation we get

`(T_p^(n+1)-T_p^n)/(Deltat)=alpha*((T_L-2*T_p+T_R)/(Deltax^2) +(T_T-2*T_p+T_B)/(Deltay^ 2))^n`

There are two methods to solve this equation

1) Explicit:- In the explicit method, there would be only one unknown. In this case it is `(T_p^(n+1)'. All other variables are at time step n which have been calculated. 

Solving for `(T_p^(n+1)', we get

`T_p^(n+1)=T_p^(n)+ ((alpha*Deltat)/(Deltax^2))*(T_L^n-2*T_p^n+T_R^n)+ ((alpha*Deltat)/(Deltay^2))*(T_T^n-2*T_p^n+T_B^n)`

In this equation we have the one unknown  `(T_p^(n+1)' on the LHS. All the parameters on RHS are at time step n , a known quantity.

In the above equation we can denote as follows.

' k1=((alpha*Deltat)/(Deltax^2))'

' k2=((alpha*Deltat)/(Deltay^2))'

'k1' and 'k2' are called the CFL numbers for the equation

For the numerical explicit diffusion equation  to be stable, it must satisfy the below criteria

`k1+k2<= 0.5`

In order to solve the transient explicit equation, we solve for Temperature values at each node using values from previous time step. We then run the loop for sufficiently large time steps to obtain the solution. No iterative solver is required.

Code for transient equation

clear all
close all
clc

nx=50;
ny=50;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
alpha=4e-3;
dt=0.01;
k1=(alpha*dt)/(dx)^2;
k2=(alpha*dt)/(dy)^2;
T=ones(nx,ny);
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
Told=T;
T(1,1)=(400+900)/2;
T(end,1)=(400+600)/2;
T(1,end)=(900+800)/2;
T(end,end)=(800+600)/2;
% k=2*((dx^2+dy^2)/(dx^2*dy^2));

Tol=1e-4;
error1=9e9;


for nt=1:15000
    tic;
    for i=2:nx-1
        for j=2:ny-1
            m=k1*(Told(i+1,j)+Told(i-1,j)-2*Told(i,j));
            n=k2*(Told(i,j+1)+Told(i,j-1)-2*Told(i,j));
            T(i,j)=m+n+Told(i,j) ;
        end
    end
    time1=toc;
Told=T;
end
figure(1);
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit1=sprintf('Simulation Time=%d s',time1)
title({'Solution for transient explicit '})
colorbar
pause(0.3)

 

Plot for temperature profile

Parameters considered

No of points along x=50

No of points along y=50

dx=dy=0.0204

thermal diffusivity 'alpha=4e-3'

dt=0.01

CFL along x 'k1=0.0960'

CFL along y 'k2=0.0960'

'k1+k2=0.192' which is less than 0.5. We expect a solution

No of time steps= 150000

Plot

 

Plot when 'k1+k2 >0.5'. The solution is unstable

Implicit Method

In this method, we have more than one unknown element

`(T_p^(n+1)-T_p^n)/(Deltat)=alpha*((T_L-2*T_p+T_R)/(Deltax^2) +(T_T-2*T_p+T_B)/(Deltay^ 2))^(n+1)`

`T_p^(n+1)=T_p^(n)+ ((alpha*Deltat)/(Deltax^2))*(T_L^(n+1)-2*T_p^(n+1)+T_R^(n+1))+ ((alpha*Deltat)/(Deltay^2))*(T_T^(n+1)-2*T_p^(n+1)+T_B^(n+1))`

In the above equation we can denote as follows.

' k1=((alpha*Deltat)/(Deltax^2))'

' k2=((alpha*Deltat)/(Deltay^2))'

Collecting terms related to `T_p^(n+1)' on LHS we get

`'T_p^(n+1)*[1+2k1+2k2]= T_p^n+(k1*(T_L^(n+1)+T_R^(n+1)) +k2*(T_T^(n+1)+T_B^(n+1)) )`

Solving for `'T_p^(n+1)'

`'T_p^(n+1)*= [1+2k1+2k2]^-1*(T_p^n+k1*(T_L^(n+1)+T_R^(n+1)) +k2*(T_T^(n+1)+T_B^(n+1)) )`

For ease of coding the terms of the above equation are labelled as

'[1+2k1+2k2]^-1= Term1'

'k1*[1+2k1+2k2]^-1= k1*Term1=Term2'

'k2*[1+2k1+2k2]^-1= k2*Term1=Term3'

The above equation is soved using the 3 iterative solvers that were used previously

1. Jacobi solver

clear all
close all
clc

nx=50;
ny=50;
x=linspace(0,1,nx);%Grid Spacing along X
y=linspace(0,1,ny);%Grid Spacing along Y

dx=x(2)-x(1);
dy=y(2)-y(1);
alpha=1.5;%Thermal Diffusivity
dt=1e-2;%Time Step
k1=(alpha*dt)/(dx)^2;%CFL along x
k2=(alpha*dt)/(dy)^2;%CFL along y
term1=(1/(1+2*k1+2*k2));
term2=k1*term1;
term3=k2*term1;
T=ones(nx,ny);
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;

 T(1,1)=(400+600)/2;
 T(end,1)=(400+800)/2;
 T(1,end)=(900+800)/2;
 T(end,end)=(800+600)/2;
 Told=T;
Tinitial=T;
Tprev=T;
 
Tol=1e-4;

jacobi_iter=1;
a=1
for nt=1:1400
    error1=9e9;
    tic;
    while(error1>Tol)
        for i=2:nx-1
            for j=2:ny-1
                H=term2*(Told(i-1,j)+Told(i+1,j));
                V=term3*(Told(i,j+1)+Told(i,j-1));
                T(i,j)=term1*Tprev(i,j)+ H+V;
            end
        end
        error1=max(max(abs(Told-T)));
        Error1(a)=error1;
        Told=T;%Spatial Temp values updated within the space iteration
        jacobi_iter=jacobi_iter+1;
        JI(a)=jacobi_iter;
        a=a+1;
            time1=toc;
    end
    Tprev=T;%After convergence the T values updated

    
end
figure(1)
colormap(jet);
contourf(x,y,T,'showText','on')
 set(gca,'YDir','reverse')
tit1=sprintf('Number of iteration=%d\n Simulation Time=%d s',jacobi_iter,time1)
title({'Jacobi Iter solver for implicit',tit1})
colorbar
pause(0.3)

Plot for Temperature 

Parameters considered

No of points along x=50

No of points along y=50

dx=dy=0.0204

thermal diffusivity 'alpha=1,5'

dt=0.01

CFL along x 'k1=0.9604'

CFL along y 'k2=0.9604'

No of Time steps = 1400

 

No of iterations to converge 24199, Simulation time 3.51e-05s

 

- Gauss seidel Method

 Code ( To be read in continuation of Implicit)

%Gauss Seidel Implicit

T=Tinitial;
Told=Tinitial;
Tprev=Tinitial;
gauss_seidel=1;
b=1;
for nt=1:1400
    error2=9e9;
    tic;
    while(error2>Tol)
        for i=2:nx-1
            for j=2:ny-1
                H=term2*(T(i-1,j)+T(i+1,j));
                V=term3*(T(i,j+1)+T(i,j-1));
                T(i,j)=term1*Tprev(i,j)+ H+V;
            end
        end
        error2=max(max(abs(Told-T)));
        Error2(b)=error2;
    Told=T;
    GS(b)=gauss_seidel;
     gauss_seidel=gauss_seidel+1;
     
     b=b+1;
    
    end
    Tprev=T;
     time2=toc;
    
   
end
figure(2)
colormap(jet);

contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse') 
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',gauss_seidel,time2)
title({'Gauss Seidel Iter solver for implicit',tit2})
colorbar
pause(0.3)

Plot for Temperature 

Parameters considered

No of points along x=50

No of points along y=50

dx=dy=0.0204

thermal diffusivity 'alpha=1,5'

dt=0.01

CFL along x 'k1=0.9604'

CFL along y 'k2=0.9604'

No of Time steps = 1400

No of Iteration to converge 14410, Simulation Time 2.7e-05 sec

-SOR

Code for SOR (To be read in continuation with Gauss Seide Implicit)

%SOR Implicit
Told=Tinitial;
Tprev=Tinitial;
T=Tinitial;
omega=1.5
SOR=1;
c=1;
for nt=1:1400
    error3=9e9;
    tic;
    while(error3>Tol)
        for i=2:nx-1
            for j=2:ny-1
                H=term2*(T(i-1,j)+T(i+1,j));
                V=term3*(T(i,j+1)+T(i,j-1));
                T(i,j)=(term1*Tprev(i,j)+ H+V)*omega+(1-omega)*Told(i,j);
            end
        end
        error3=max(max(abs(Told-T)));
        Error3(c)=error3;
    Told=T;
    SR(c)=SOR;
     SOR=SOR+1;
  c=c+1;
      time3=toc;
    end
    Tprev=T;
   
   
end
 figure(3)
colormap(jet);

contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',SOR,time3)
title({'SOR Iter solver for implicit',tit2})
colorbar

pause(0.3)


figure(4)
semilogy(JI,Error1,GS,Error2,'r',SR, Error3,'k');

legend('Jacobi','Gauss Seidel','SOR')

Plot for Temperature 

Parameters considered

No of points along x=50

No of points along y=50

dx=dy=0.0204

thermal diffusivity 'alpha=1,5'

dt=0.01

relaxation factor omega=1.5

CFL along x 'k1=0.9604'

CFL along y 'k2=0.9604'

No of Time steps = 1400

No of Iteration to converge 6614, Simulation Time 3.3e-05 sec

 

 

-Comparison of decrease of error with iteration for 3 Methods

 

 

From the above graph we note that in SOR method the rate of decrease of error is much higher compared to Gauss seidel and Jacobi. Jacobi takes much longer time since we continue to use the old guess values instead of the newly calculated better guess values.

 

%Full COde for Steady State
clear all
close all
clc

%Input Parameter
nx=50;
ny=50;
x=linspace(0,1,nx);%Grid along X
y=linspace(0,1,ny);%Grid along y
dx=x(2)-x(1);
dy=y(2)-y(1);
T=ones(nx,ny);

%Boundary Conditions
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
Tinitial=T;
Told=T;
k=2*((dx^2+dy^2)/(dx^2*dy^2));

Tol=1e-4
error1=9e9
error2=9e9
error3=9e9
jacobi_iter=1;

tic
a=1;
while (error1>Tol)
    for i=2:nx-1
        for j=2:ny-1
            m=(Told(i+1,j)+Told(i-1,j))/(dx*dx);
            n=(Told(i,j+1)+Told(i,j-1))/(dy*dy);
            T(i,j)=(m+n)/k ;
        end
    end
    error1= max(max(abs(Told-T)));
    ERROR1(a)=error1; %Plotting error vs iterations
    Told=T;
    jacobi_iter=jacobi_iter+1;
    JI(a)=jacobi_iter;%Storing iterations in counter
    time1=toc;
    a=a+1;
    

end
figure(1);
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit1=sprintf('Number of iteration=%d\n Simulation Time=%d s',jacobi_iter,time1)
title({'Jacobi solver for steady state',tit1})
colorbar

pause(0.3)

%Gauss Seidel

tic;
T=Tinitial;
Told=Tinitial;
gauss_seidel=1;
b=1;
while (error2>Tol)
    for i=2:nx-1
        for j=2:ny-1
            m=(T(i+1,j)+T(i-1,j))/(dx*dx);
            n=(T(i,j+1)+T(i,j-1))/(dy*dy);
            T(i,j)=(m+n)/k ;
        end
    end
    error2= max(max(abs(Told-T)));
      Told=T;
      ERROR2(b)= error2;
       GS(b)=gauss_seidel;
    gauss_seidel=gauss_seidel+1;
   b=b+1;
    time2=toc;

end
figure(2)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',gauss_seidel,time2)
title({'Gauss Seidel solver for steady state',tit2})
colorbar
pause(0.3)

%Successive Over Relaxation
SOR=1;

T=Tinitial
Told=Tinitial;
omega=1.5;
tic

c=1;
while (error3>Tol)
    for i=2:nx-1
        for j=2:ny-1
            m=(T(i+1,j)+T(i-1,j))/(dx*dx);
            n=(T(i,j+1)+T(i,j-1))/(dy*dy);
            
            T(i,j)=(1-omega)*Told(i,j) +omega*((m+n)/k );
        end
    end
    error3= max(max(abs(Told-T)));
    ERROR3(c)= error3;
    Told=T;
    SR(c)=SOR;
    SOR=SOR+1;
    time3=toc;
    c=c+1;

end
figure(3)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit3=sprintf('Number of iteration=%d\n Simulation Time=%d s',SOR,time3)
title({'SOR solver for steady state',tit3})
colorbar
pause(0.3)

figure(4)
semilogy(JI,ERROR1,GS,ERROR2,'r',SR, ERROR3,'k');

legend('Jacobi','Gauss Seidel','SOR')
xlabel('No of iterations')
ylabel('log of absolute error')
title ('Comaparison of rate of decrease of error with iterations for three methods')

%Full Code for Transient explicit

clear all
close all
clc

nx=50;
ny=50;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
alpha=4e-3;
dt=0.01;
k1=(alpha*dt)/(dx)^2;
k2=(alpha*dt)/(dy)^2;
T=ones(nx,ny);
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
Told=T;
T(1,1)=(400+900)/2;
T(end,1)=(400+600)/2;
T(1,end)=(900+800)/2;
T(end,end)=(800+600)/2;
% k=2*((dx^2+dy^2)/(dx^2*dy^2));

Tol=1e-4;
error1=9e9;


for nt=1:15000
    tic;
    for i=2:nx-1
        for j=2:ny-1
            m=k1*(Told(i+1,j)+Told(i-1,j)-2*Told(i,j));
            n=k2*(Told(i,j+1)+Told(i,j-1)-2*Told(i,j));
            T(i,j)=m+n+Told(i,j) ;
        end
    end
    time1=toc;
Told=T;
end
figure(1);
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit1=sprintf('Simulation Time=%d s',time1)
title({'Solution for transient explicit '})
colorbar
pause(0.3)

%Full COde for Implicit

clear all
close all
clc

nx=50;
ny=50;
x=linspace(0,1,nx);%Grid Spacing along X
y=linspace(0,1,ny);%Grid Spacing along Y

dx=x(2)-x(1);
dy=y(2)-y(1);
alpha=1.5;%Thermal Diffusivity
dt=1e-2;%Time Step
k1=(alpha*dt)/(dx)^2;%CFL along x
k2=(alpha*dt)/(dy)^2;%CFL along y
term1=(1/(1+2*k1+2*k2));
term2=k1*term1;
term3=k2*term1;
T=ones(nx,ny);
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;

 T(1,1)=(400+600)/2;
 T(end,1)=(400+800)/2;
 T(1,end)=(900+800)/2;
 T(end,end)=(800+600)/2;
 Told=T;
Tinitial=T;
Tprev=T;
 
Tol=1e-4;

jacobi_iter=1;
a=1
for nt=1:1400
    error1=9e9;
    tic;
    while(error1>Tol)
        for i=2:nx-1
            for j=2:ny-1
                H=term2*(Told(i-1,j)+Told(i+1,j));
                V=term3*(Told(i,j+1)+Told(i,j-1));
                T(i,j)=term1*Tprev(i,j)+ H+V;
            end
        end
        error1=max(max(abs(Told-T)));
        Error1(a)=error1;
        Told=T;%Spatial Temp values updated within the space iteration
        jacobi_iter=jacobi_iter+1;
        JI(a)=jacobi_iter;
        a=a+1;
            time1=toc;
    end
    Tprev=T;%After convergence the T values updated

    
end
figure(1)
colormap(jet);
contourf(x,y,T,'showText','on')
 set(gca,'YDir','reverse')
tit1=sprintf('Number of iteration=%d\n Simulation Time=%d s',jacobi_iter,time1)
title({'Jacobi Iter solver for implicit',tit1})
colorbar
pause(0.3)

%Gauss Seidel Implicit

T=Tinitial;
Told=Tinitial;
Tprev=Tinitial;
gauss_seidel=1;
b=1;
for nt=1:1400
    error2=9e9;
    tic;
    while(error2>Tol)
        for i=2:nx-1
            for j=2:ny-1
                H=term2*(T(i-1,j)+T(i+1,j));
                V=term3*(T(i,j+1)+T(i,j-1));
                T(i,j)=term1*Tprev(i,j)+ H+V;
            end
        end
        error2=max(max(abs(Told-T)));
        Error2(b)=error2;
    Told=T;
    GS(b)=gauss_seidel;
     gauss_seidel=gauss_seidel+1;
     
     b=b+1;
    
    end
    Tprev=T;
     time2=toc;
    
   
end
figure(2)
colormap(jet);

contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse') 
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',gauss_seidel,time2)
title({'Gauss Seidel Iter solver for implicit',tit2})
colorbar
pause(0.3)

%SOR Implicit
Told=Tinitial;
Tprev=Tinitial;
T=Tinitial;
omega=1.5
SOR=1;
c=1;
for nt=1:1400
    error3=9e9;
    tic;
    while(error3>Tol)
        for i=2:nx-1
            for j=2:ny-1
                H=term2*(T(i-1,j)+T(i+1,j));
                V=term3*(T(i,j+1)+T(i,j-1));
                T(i,j)=(term1*Tprev(i,j)+ H+V)*omega+(1-omega)*Told(i,j);
            end
        end
        error3=max(max(abs(Told-T)));
        Error3(c)=error3;
    Told=T;
    SR(c)=SOR;
     SOR=SOR+1;
  c=c+1;
      time3=toc;
    end
    Tprev=T;
   
   
end
 figure(3)
colormap(jet);

contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',SOR,time3)
title({'SOR Iter solver for implicit',tit2})
colorbar

pause(0.3)


figure(4)
semilogy(JI,Error1,GS,Error2,'r',SR, Error3,'k');
xlabel('ITeration')
ylabel('log of absolute error')
title('Comparison of error with iterations')
legend('Jacobi','Gauss Seidel','SOR')