Eliptic Partial Differential Equations 
Direct-Nicolson Method

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 direct method 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

l=2.4;
w=3.0;
TL=75.0;
TR=100.0;
TT=300.0;
TB=50.0;
dx=0.6;
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];

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];

Explicit Method

Defining the coefficient  matrix

A=Table[0,{i,1,(m-1)*(n-1)},{j,1,(m-1)*(n-1)}];

Defining the right hand side matrix

B=Table[0,{i,1,(m-1)*(n-1)},{j,1,1}];

Loop to create the left hand side coefficient matrix

count=0;
For[i=1,i≤(m-1)*(n-1),i=i+1,
j=i;
count=count+1;
A[[i,j]]=-4;
If[count≠1,A[[i,j-1]]=1];
If[count≠n-1,A[[i,j+1]]=1];
If[j+n-1≤(m-1)*(n-1),A[[i,j+n-1]]=1];
If[j-n+1>0,A[[i,j-n+1]]=1];
If[count==n-1,count=0];
]

Loop to create the right hand side matrix

count=0;
For[i=1,i≤m-1,i=i+1,
For[j=1,j≤n-1,j=j+1,
count=count+1;
temp1=0;
temp2=0;
If[i-1==0,temp1=TL, If[i+1==m,temp1=TR]  ];
If[j-1==0,temp2=TT, If[j+1==n,temp2=TB]  ];
B[[count,1]]=-(temp1+temp2);
]
]

Solving for temperatures of the internal nodes

sol=Table[Inverse[A].B];

Creating a matrix with temperatures for internal nodes

Tinternal=Table[0,{i,1,m-1},{j,1,n-1}];
count=0;
For[i=1,i≤m-1,i=i+1,
For[j=1,j≤n-1,j=j+1,
count=count+1;
Tinternal[[i,j]]=sol[[count,1]];
]
]
Tinternal=Transpose[Tinternal];

Creating a matrix,T to store the temperatures of the internal and boundary nodes

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

Assigning the boundary temperatures to the temperature matrix,T

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

Assigning the internal node temperatures to the temperature matrix,T

For[i=1,i≤n-1,i=i+1,
For[j=1,j≤m-1,j=j+1,
T[[i+1,j+1]]=Tinternal[[i,j]];
]
]

Displaying the temperature matrix, T

Print["Temperature profile, T="];Print[TableForm[Tdisplay]];
MatrixPlot[T,ColorFunction->"Rainbow"]