计算科学 - 有限差异与元素：准确性和实施 - 吾爱随笔录

我正在尝试以数值方式求解 2D Poisson 方程：

$\frac{\partial ^2 \phi}{\partial x^2} + \frac{\partial ^2 \phi}{\partial y^2} = 1$

与狄利克雷边界条件 $\phi = 0$ .

我在 6 × 6 点的规则网格上使用了有限元和有限差分方法。模型尺寸为 1 × 1 m。两种方法的结果似乎完全不同（我还测试了两种方法之间具有相似差异的较大网格）。

结果有限元： $\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.        ]]

结果有限差分： $\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.        ]]

此外，如果我在左侧设置无通量边界条件，则网格第一列和第二列的结果为：

结果有限元： $\phi =$

[[ 0.          0.        ]
 [-0.07453021 -0.0730567 ]
 [-0.11114969 -0.10876554]
 [-0.11114969 -0.10876554]
 [-0.07453021 -0.0730567 ]
 [ 0.          0.        ]]

结果有限差分： $\phi =$

[[ 0.          0.        ]
 [-0.06977283 -0.06977283]
 [-0.1034833  -0.1034833 ]
 [-0.1034833  -0.1034833 ]
 [-0.06977283 -0.06977283]
 [ 0.          0.        ]]

有限差分法的结果表明，在左边界处实际上存在零通量。然而，这不是有限元方法的情况。

从这些结果中，我得到的印象是有限差分更准确。这是正确的（一般）？还是我错误地实现了数值方法（参见下面的代码）？由于这两种方法都基于不同的假设，因此我得出了不同的结果。然而，这些结果似乎差异太大而不能成立。

在这里，我提供了我实现的 Python 代码，用于使用有限元和有限差分来求解泊松方程。

##################################################################
### 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