# Finite differences vs. elements: Accuracy and implementation

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

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

-
Why not try the Method of Manufactured solutions and compared the error when you know the exact solution? There's no way to tell from your answers which is closer to the true solution. – Bill Barth Jun 7 '14 at 16:05

In the finite element method, you impose the no-flux boundary conditions only weakly. In other words, the flux you get (say, $g=\partial u_h/\partial n$) is not zero at every point. All you can say is that $(g,\varphi_i)_{\partial\Omega}=0$ for all test functions. In the limit $h\rightarrow 0$, this guarantees that the flux is zero, but not on any given mesh -- there, it is simply zero with respect to some test functions. – Wolfgang Bangerth Jun 8 '14 at 0:27