4
$\begingroup$

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
Improve this question
$\endgroup$
1
  • 3
    $\begingroup$ 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. $\endgroup$
    – Bill Barth
    Jun 7 '14 at 16:05
6
$\begingroup$

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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$
    – jan
    Jun 7 '14 at 10:46
  • $\begingroup$ 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? $\endgroup$
    – Godric Seer
    Jun 7 '14 at 13:52
  • $\begingroup$ 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. $\endgroup$
    – jan
    Jun 7 '14 at 16:09
  • $\begingroup$ 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. $\endgroup$
    – Godric Seer
    Jun 7 '14 at 17:02
  • 3
    $\begingroup$ 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. $\endgroup$
    – Wolfgang Bangerth
    Jun 8 '14 at 0:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.