# 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

# 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

• 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. Jun 7, 2014 at 16:05

For a 6x6 grid, those are about the error differences I would expect from two different methods. You have to realize that a 6x6 grid is a very coarse grid, even for a simple problem like yours. As long as you see the two solutions converge towards each other as you refine your grid, there is likely no implementation error. Finite-difference has no general advantage over finite-element, although in some problems it may be more accurate up to a constant factor (and the opposite is true for other problems).

As far as the Neumann condition, if you formulated it correctly in your FE case, then you may not see column 1 equal to column 2. What you should actually check is that the approximating functions defined over your elements have zero derivative at the boundary. However, if you are using first order elements (which I am assuming you are), then your function are linear and my intuition would be that the two columns would be equal. You will just have to check exactly how you implemented your boundary condition.

• Thanks for your reply. I increased the number of nodes and the error differences between both methods decreases. I'm still confused about the no-flux boundary condition. Since this kind of boundary condition become part of the weak formulation of the problem, I did not change the equation for elemets of the left side of the grid (this is where there is the no-flux boundary). Only in case of the dirichlet boundaries, the equations were adjusted. Would this be the right implementation?
– jan
Jun 7, 2014 at 10:46
• Finite element is by far my weakest method, so take this with a grain of salt. If I remember correctly, on areas where you have tractions defined (derivatives, not values), that ends up on the right hand side of your system. So I think you are right, that you would just not specify the value at those nodes, and leave those nodes to take whatever value they need. Did the two columns become closer when you refined the mes? Jun 7, 2014 at 13:52
• Yes, both columns become closer after refining the grid. Based upon the results from the finite difference method, I was just expecting that the result (the two columns) were similar close in value. I guess I'll have to think about the assumptions of both methods to figure why it is not as I expected.
– jan
Jun 7, 2014 at 16:09
• The important thing to remember is that with finite element, you aren't using your standard finite difference to define that boundary condition. In reality you are defining functions over the elements, not just the values at points. You need to work with these functions, not the points. Jun 7, 2014 at 17:02
• 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. Jun 8, 2014 at 0:27