Symbolic generating of system equations for 2D regular grid for solving Laplace equation using finite difference method
By Nasser M. Abbasi
updated oct 28,2010

Introduction

When solving (∂_(x,x))u+(∂_(y,y))u=f(x,y) using finite difference method, in order to make it easy to see the internal structure of the A matrix using the standard 5 points Laplacian scheme, the following is a small function which generates the symbolic format of these equations for  a given N, the number of grid points on one edge.  At the end of this note, the system equations are generated for N=4,5,6,7,8. One can see the form of the A matrix with the dominant diagonal and the corresponding bands. It is mostly a sparse matrix.

The indexing method used is that described in the class.

Define U at each grid point

n=4;
u=Array[(Subscript[U, ##])&,{n,n},{0,0}];
Style[MatrixForm@Reverse@u,18]

Define the force vector

f=Array[(Subscript[F, ##])&,{n,n},{0,0}];
Style[MatrixForm@Reverse@f,18]

Needs["Notation`"]
Notation[(NotationTemplateTag[(Subscript[f, i_,j_])]) ⟺ (NotationTemplateTag[f[[i_+1,j_+1]]])]
Notation[(NotationTemplateTag[(Subscript[u, i_,j_])]) ⟺ (NotationTemplateTag[u[[i_+1,j_+1]]])]

eq[u_,i_,j_,h_]:=Module[{},
(((Subscript[u, i-1,j])+(Subscript[u, i+1,j])+(Subscript[u, i,j-1])+(Subscript[u, i,j+1])-4 (Subscript[u, i,j]))/(h^2))
];

Draw the grid with the unknown above at each point

nRow=4;nCol=4;
makeGrid[nRow_,nCol_,u_]:=Module[{i,j},
internalGrid=Table[{
Line[{{i-1,j},{i+1,j}}],
Line[{{i,j-1},{i,j+1}}],
{PointSize[Large],Point[{i,j}]},
Text[(Subscript[u, i,j]),{i,j},{1.5,1.5}]
},
{j,1,nCol-2},{i,1,nRow-2}
];
boundary={
Line[{{0,0},{nCol-1,0}}],
Line[{{0,0},{0,nRow-1}}],
Line[{{nCol-1,nRow-1},{nCol-1,0}}],
Line[{{nCol-1,nRow-1},{0,nRow-1}}]};
Graphics[boundary~Join~internalGrid]
]
makeGrid[nRow,nCol,u]

Generate the discrete equations at each of the internal grid points

n=4;
Style[MatrixForm[eqs=Flatten@Table[eq[u,i,j,h]==(Subscript[f, i,j]),{j,n-2},{i,n-2}]],22]

List the unknowns

Style[MatrixForm[unknowns=Flatten@Table[(Subscript[u, i,j]),{j,n-2},{i,n-2}]],22]

Generate the equations of the form AU=F

unknowns=Flatten@Table[(Subscript[u, i,j]),{j,n-2},{i,n-2}];
{b,A}=Normal@CoefficientArrays[eqs,unknowns];
Style[Grid[{{(1/(h^2)),MatrixForm[(h^2)A],MatrixForm@unknowns,"=",MatrixForm@b}}],22]

For homegenous Boundary conditions

Table[(Subscript[u, 0,j])->0,{j,0,3}]~Join~Table[(Subscript[u, i,0])->0,{i,0,3}]~Join~Table[(Subscript[u, 3,j])->0,{j,0,3}]~Join~Table[(Subscript[u, i,3])->0,{i,0,3}];
{b0,A}=Normal@CoefficientArrays[eqs,unknowns]/.%;
Style[Grid[{{(1/(h^2)),MatrixForm[(h^2)A],MatrixForm@unknowns,"=",MatrixForm@b0}}],22]

set the boundary conditions. Assume left side is U=α,Right side U=β, bottom side U=γ, top side U=η, then the above becomes

b=b/.Subscript[U,0,any_]->α; (*left*)
b=b/.Subscript[U,any_,3]->β; (*right*)
b=b/.Subscript[U,any_,0]->γ; (*bottom*)
b=b/.Subscript[U,3,any_]->η;(*top*)

Display the equations again

Style[Grid[{{(1/(h^2)),MatrixForm[(h^2)A],MatrixForm@unknowns,"=",MatrixForm@b}}],22]

Now the system can be solved for the unknowns U , given the force F values.

Below is the system equations generated for N=3,4,5,6,7,8. Put the above code into one function to use it all the time

generateSystemEquations[n_,U_,F_,h_,α_,β_,γ_,η_]:=Module[{m,f,i,j,eqs,unknowns,b,A,u},
u=Array[(Subscript[U, ##])&,{n,n},{0,0}];
f=Array[(Subscript[F, ##])&,{n,n},{0,0}];
eqs=Flatten@Table[eq[u,i,j,h]==(Subscript[f, i,j]),{j,n-2},{i,n-2}];
unknowns=Flatten@Table[(Subscript[u, i,j]),{j,n-2},{i,n-2}];
{b,A}=Normal@CoefficientArrays[eqs,unknowns];
b=b/.Subscript[U,0,any_]->α; (*left*)
b=b/.Subscript[U,any_,n-1]->β; (*right*)
b=b/.Subscript[U,any_,0]->γ; (*bottom*)
b=b/.Subscript[U,n-1,any_]->η;(*top*)
Print[makeGrid[n,n,u]];
Print[Style[Row[{(1/(h^2)),MatrixForm[(h^2)A],MatrixForm@unknowns,"=",MatrixForm@b},Alignment->Left,ImageSize->Full],18]]
]

Clear[U,F,h,α,β,γ,η];
Do[{generateSystemEquations[i,U,F,h,α,β,γ,η],Print@Graphics[{Thick,Line[{{0,0},{10,0}}]}]},{i,4,6}]

For homogenous boundary conditions

generateSystemEquations[n_,U_,F_,h_,α_,β_,γ_,η_]:=Module[{m,f,i,j,eqs,unknowns,b,A,u},
u=Array[(Subscript[U, ##])&,{n,n},{0,0}];
f=Array[(Subscript[F, ##])&,{n,n},{0,0}];
eqs=Flatten@Table[eq[u,i,j,h]==(Subscript[f, i,j]),{j,n-2},{i,n-2}];
unknowns=Flatten@Table[(Subscript[u, i,j]),{j,n-2},{i,n-2}];
{b,A}=Normal@CoefficientArrays[eqs,unknowns];
b=b/.Subscript[U,0,any_]->0; (*left*)
b=b/.Subscript[U,any_,n-1]->0; (*right*)
b=b/.Subscript[U,any_,0]->0; (*bottom*)
b=b/.Subscript[U,n-1,any_]->0;(*top*)
Print[makeGrid[n,n,u]];
Print[Style[Row[{(1/(h^2)),MatrixForm[(h^2)A],MatrixForm@unknowns,"=",MatrixForm@b},Alignment->Left,ImageSize->Full],18]]
]
Do[{generateSystemEquations[i,U,F,h,α,β,γ,η],Print@Graphics[{Thick,Line[{{0,0},{10,0}}]}]},{i,4,6}]