I am trying to solve the 2D Poisson equation numerically:
$ \frac{\partial ^2 \phi}{\partial x^2} + \frac{\partial ^2 \phi}{\partial y^2} = 1 $
with the Dirichlet boundary condition $\phi = 0$.
I used the finite element and finite difference approach on a regular grid of 6 × 6 points. The model size is 1 × 1 m. The results appear to be quite different of both methods (I also tested larger grids with similar differences between both methods).
Results finite element: $\phi = $
[[ 0. 0. 0. 0. 0. 0. ]
[ 0. -0.03631579 -0.04973684 -0.04973684 -0.03631579 0. ]
[ 0. -0.04973684 -0.07105263 -0.07105263 -0.04973684 0. ]
[ 0. -0.04973684 -0.07105263 -0.07105263 -0.04973684 0. ]
[ 0. -0.03631579 -0.04973684 -0.04973684 -0.03631579 0. ]
[ 0. 0. 0. 0. 0. 0. ]]
Results finite difference: $\phi = $
[[ 0. 0. 0. 0. 0. 0. ]
[ 0. -0.03333333 -0.04666667 -0.04666667 -0.03333333 0. ]
[ 0. -0.04666667 -0.06666667 -0.06666667 -0.04666667 0. ]
[ 0. -0.04666667 -0.06666667 -0.06666667 -0.04666667 0. ]
[ 0. -0.03333333 -0.04666667 -0.04666667 -0.03333333 0. ]
[ 0. 0. 0. 0. 0. 0. ]]
Furtheremore, if I set a no-flux boundary condition on the left side, the results for the first and second column of the grid are:
Results finite element: $\phi = $
[[ 0. 0. ]
[-0.07453021 -0.0730567 ]
[-0.11114969 -0.10876554]
[-0.11114969 -0.10876554]
[-0.07453021 -0.0730567 ]
[ 0. 0. ]]
Results finite difference: $\phi = $
[[ 0. 0. ]
[-0.06977283 -0.06977283]
[-0.1034833 -0.1034833 ]
[-0.1034833 -0.1034833 ]
[-0.06977283 -0.06977283]
[ 0. 0. ]]
The results of the finite difference method show there is actually a zero flux at the left boundary. However, this is not the case for the finite element approach.
From these results, I get the impression that the finite difference is more accurate. Is this correct (in general)? Or have I implemented the numerical approach incorrectly (see code below)? Since both methods are based on different assumptions, I expeted different results. However, these results seem to be too different to be true.
Here, I provide the Python code I implemented to solve the Poisson equation using finite elements and finite differences.
##################################################################
### IMPORT ###
##################################################################
from numpy import zeros,sqrt,dot,transpose,sqrt
from numpy.linalg import det,inv
from scipy.sparse.linalg import spsolve
from scipy.sparse import csc_matrix
##################################################################
### SETUP ###
##################################################################
nnx = 6 # number of nodes - x axis
nny = 6 # number of nodes - y axis
np = nnx*nny # total number of nodes
nelx = nnx-1 # number of elements - x axis
nely = nny-1 # number of elements - y axis
nel = nelx*nely # total number of elements
Lx = 1.0 # x axis goes from 0 to Lx
Ly = 1.0 # x axis goes from 0 to Lx
xstp = Lx/(nnx-1) # x step size
ystp = Ly/(nnx-1) # y step size
x=zeros((np,1))
y=zeros((np,1))
ind=-1
for j in range(nny):
for i in range(nnx):
ind=ind+1
x[ind,0]=i*Lx/nelx
y[ind,0]=j*Ly/nely
##################################################################
#****************************************************************#
#FINITE ELEMENT APPROACH
#****************************************************************#
##################################################################
##################################################################
### CONNECTIVITY OF NODES FOR EACH ELEMENT ###
##################################################################
icon=zeros((4,nel))
ind=-1
eind=-1
for j in range(nny):
for i in range(nnx):
ind=ind+1
if j==nny-1 or i==nnx-1:
continue
eind += 1
icon[0,eind]=ind
icon[1,eind]=ind+1
icon[2,eind]=ind+1+nnx
icon[3,eind]=ind+nnx
##################################################################
### BOUNDARY CONDITIONS SETUP ####
##################################################################
bc_fix=zeros((np,1))
bc_val=zeros((np,1))
for i in range(np):
if x[i,0]==0.0:
bc_fix[i,0] = 1
bc_val[i,0] = 0.0
if y[i,0]==0.0:
bc_fix[i,0] = 1
bc_val[i,0] = 0.0
if x[i,0]==Lx:
bc_fix[i,0] = 1
bc_val[i,0] = 0.0
if y[i,0]==Ly:
bc_fix[i,0] = 1
bc_val[i,0] = 0.0
##################################################################
### ASSEMBLY ###
##################################################################
A = zeros((np,np)) # GLOBAL MATRIX - LHS
B = zeros((np,1)) # GLOBAL MATRIX - RHS
# Weights for quadtratic integrations
wgts = [1.0]*4
# Integration points
intpt_x = [-1.0/sqrt(3),-1.0/sqrt(3), 1.0/sqrt(3), 1.0/sqrt(3)]
intpt_y = [-1.0/sqrt(3), 1.0/sqrt(3),-1.0/sqrt(3), 1.0/sqrt(3)]
for iel in range(nel): # loop over each element
Ael=zeros((4,4)) # element matrix
Bel=zeros((4,1)) # element matrix
for i in range(4): # loop over each integration point
wq=wgts[i]
rq=intpt_x[i]
sq=intpt_y[i]
# Shape Function
N = zeros((4,1))
N[0,0]=0.25*(1.0-rq)*(1.0-sq)
N[1,0]=0.25*(1.0+rq)*(1.0-sq)
N[2,0]=0.25*(1.0+rq)*(1.0+sq)
N[3,0]=0.25*(1.0-rq)*(1.0+sq)
# Shape function derivatives
dNdrs = zeros((4,2))
dNdrs[0,0] = - 0.25*(1.0-sq)
dNdrs[1,0] = + 0.25*(1.0-sq)
dNdrs[2,0] = + 0.25*(1.0+sq)
dNdrs[3,0] = - 0.25*(1.0+sq)
dNdrs[0,1] = - 0.25*(1.0-rq)
dNdrs[1,1] = - 0.25*(1.0+rq)
dNdrs[2,1] = + 0.25*(1.0+rq)
dNdrs[3,1] = + 0.25*(1.0-rq)
# Calculate Jacobian
cord = zeros((2,4)) # cordinates of element
for j in range(4):
cord[0,j] = x[icon[j,iel]]
cord[1,j] = y[icon[j,iel]]
J = dot(cord,dNdrs) # jacobian
detJ = det(J) # determinant
invJ = inv(J) # inverse jacobian
# Local Derivatives
dNdrs_l = dot(dNdrs,invJ)
# Create Element Matrix
Ael -= dot(dNdrs_l,transpose(dNdrs_l))*detJ*wq
Bel += N*detJ*wq
# Update Global Matrix
for k1 in range(4):
ik1=icon[k1,iel]
for k2 in range(4):
ik2=icon[k2,iel]
A[ik1,ik2] += Ael[k1,k2]
B[ik1]=B[ik1]+Bel[k1]
##################################################################
# SET BOUNDARY CONDITIONS
##################################################################
for i in range(np):
if bc_fix[i] == 1:
for j in range(np):
B[j]=B[j]-A[i,j]*bc_val[i]
A[i,j]=0.0
A[j,i]=0.0
A[i,i]=1.0
B[i]=bc_val[i]
##################################################################
# SOLVE ...
##################################################################
A = csc_matrix(A)
S_fe = spsolve(A,B)
S_fe = S_fe.reshape(nny,nnx)
print S_fe
##################################################################
#****************************************************************#
#FINITE DIFFERENCE APPROACH
#****************************************************************#
##################################################################
A = zeros((np,np))
B = zeros((np,1))
k = -1
for i in range(nny):
for j in range(nnx):
k += 1
if i==0 or i==nny-1 or j==0 or j==nnx-1:
A[k,k] = 1.0
B[k,0] = 0.0
else:
A[k,k-nny] = 1.0/xstp**2
A[k,k-1 ] = 1.0/ystp**2
A[k,k ] = -2.0/xstp**2 - 2.0/ystp**2
A[k,k+1 ] = 1.0/ystp**2
A[k,k+nny] = 1.0/xstp**2
B[k,0 ] = 1.0
##################################################################
# SOLVE ...
##################################################################
A = csc_matrix(A)
S_fd = spsolve(A,B)
S_fd = S_fd.reshape(nnx,nny)
print S_fd