Eliptic Partial Differential Equations 
Gauss-Siedel Method with Over Relaxation

Sri Harsha Garapati, Autar Kaw
University of South Florida
United States of America
kaw@eng.usf.edu

Introduction

This worksheet demonstrates the use of Mathematica to illustrate the Gauss-Siedel method with over relaxation of solving eliptic partial differential equations.

Problem Statement

Mathematica worksheet to solve for the temperature distribution of a rectangular plate with a length of 2.4 m and width of 3.0 m. The plate is subjected to dirichlet temperature boundary conditions on the foursides of the plate

Input

The following are the inputs to this problem
      a) Length of the plate, l
      b) Width of the plate, w
      c) Temperature on the left end of the plate, TL
      d) Temperature on the right end of the plate, TR
      e) Temperature on the top end of the plate, TT
      f) Temperature on the bottom end of the plate, TB
      g) grid length of the square mesh, dx
      h) Allowable tolerance for error in percentage
      i) Over relaxation factor

l=2.4;
w=3.0;
TL=75.0;
TR=100.0;
TT=300.0;
TB=50.0;
dx=0.6;
tol=1.0;
lambda=1.4;
Print["Length of the plate=",l];
Print["width of the plate=",w];
Print["Temperature on the left end of the plate=",TL];
Print["Temperature on the right end of the plate=",TR];
Print["Temperature on the top end of the plate=",TT];
Print["Temperature on the bottom end of the plate=",TB];
Print["Allowable Tolerance=",tol];
Print[" Over relaxation factor=",lambda];

Calculations

a) Square mesh so dx=dy
b) Number of elements along the length of the plate, m
c) Number of elements along the width of the plate, n

dy=dx;
m=l/dx;
n=w/dy;
Print["Number divisions along the length of the plate=",m];
Print["Number of divisions along the width of the plate=",n];

Gauss-Siedel Method

Defining temperature  matrix ,T

T=Table[0,{i,1,n+1},{j,1,m+1}];
epsa=Table[0,{i,1,n-1},{j,1,m-1}];
Tdisplay=Table[0,{i,1,n+1},{j,1,m+1}];
epsadisplay=Table[0,{i,1,n-1},{j,1,m-1}];

Loop assigning the boundary conditions to the temperature matrix

For[i=1,i≤n+1,i=i+1,
For[j=1,j≤m+1,j=j+1,
If[i==1,T[[i,j]]=TB  ];
If[i==n+1,T[[i,j]]=TT  ];
If[j==1,T[[i,j]]=TL  ];
If[j==m+1,T[[i,j]]=TR  ];
]
]

Loop converting the temperature matrix to display form

For[i=1,i≤n+1,i=i+1,
For[j=1,j≤m+1,j=j+1,
Tdisplay[[i,j]]=T[[n+1-i+1,j]];
]
]
iteration=0;
Print["Iteration Number=",iteration];
Print["Temperature profile, T="];Print[TableForm[Tdisplay]];

programming Gauss-Siedel method

count=(m-1)*(n-1);
While[count>0,
count=0;
For[i=2,i≤n,i=i+1,
For[j=2,j≤m,j=j+1,
old=T[[i,j]];
T[[i,j]]=(T[[i+1,j]]+T[[i-1,j]]+T[[i,j+1]]+T[[i,j-1]])/4;
T[[i,j]]=(lambda*T[[i,j]])+((1-lambda)*old);
epsa[[i,j]]=Abs[(T[[i,j]]-old)/T[[i,j]]*100];
If[epsa[[i,j]]>tol,count=count+1];
]
]

iteration=iteration+1;
For[i=1,i≤n+1,i=i+1,
For[j=1,j≤m+1,j=j+1,
Tdisplay[[i,j]]=T[[n+1-i+1,j]];
]
]
For[i=2,i≤n,i=i+1,
For[j=2,j≤m,j=j+1,
epsadisplay[[i,j]]=epsa[[n+1-i+1,j]];
]
]
Print["Iteration Number=",iteration];
Print["Temperature profile, T="];Print[TableForm[Tdisplay]];
]