Python实现非正则三角网格离散拉普拉斯算子特征值绝对值偏大问题排查

现象说明 #

基于scipy.spatial.Delaunay生成非正则三角网格，完成网格上离散拉普拉斯算子的实现（算子计算对应代码中名为M的函数）。代码运行无报错，但计算结果存在偏差：输出的特征向量结果正确，但得到的特征值绝对值显著偏高。实现逻辑参考公开资料中离散拉普拉斯算子条目下「网格拉普拉斯」章节给出的计算公式。

配套函数逻辑说明 #

• in_curve：返回满足f(x,y,*fargs) < 0条件的曲线内部点，测试示例中为方形区域内点生成逻辑。
• triang：返回经点加密补充的三角网格采样点，调用专用边界点生成函数（测试示例中为rect_rim函数）获取区域边界点以提升计算精度。

核心实现代码 #

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.spatial

def rect_drum(L,H,U):
    vals = []
    val = 0
    k = 1
    l = 1
    while val >= -U:
        while val >= -U:
            val = -np.pi**2*((k/L)**2+(l/H)**2)
            if val >= -U:
                vals.append(val)
            l += 1
        l = 1
        k += 1
        val = -np.pi**2*((k/L)**2+(l/H)**2)
    return np.array(vals)


def count_vp(tab,U):
    # 统计数组tab中大于等于-U的特征值数量
    return tab[tab>=-U]

def in_curve(f,fargs,shape,a):
    points = [] # 曲线内部点集
    for j in range(shape[0]):
        for i in range(shape[1]):
            if f(i*a,j*a,*fargs) < 0:
                points.append([i*a,j*a])

    return np.array(points)

def triang(points,a,f,fargs,bord):
    tri_points = points.copy()
    tri_points[:,1] *= np.sqrt(3)
    tri_points2 = np.vstack((points,bord))
    tri_points2[:,1] *= np.sqrt(3)
    tri_points2[:,0] += a/2
    tri_points2[:,1] += np.sqrt(3)/2*a
    fin = np.vstack((tri_points,tri_points2))
    i = 0
    eps = 0.01
    while i < len(fin):
        if f(fin[i,0]+eps,fin[i,1]+eps,*fargs) > 0:
            fin = np.delete(fin,i,0)
            i -= 1
        i += 1
    return np.vstack((fin,bord)),len(fin),len(bord)


def tri_ang(points,ind,p0):
    # 以p0为极点按极角排序点集
    vec=np.arctan2((points-p0)[:,1],(points-p0)[:,0])
    values = []
    dtype = [('val',float),('n',int)]
    
    for i in range(len(vec)):
        values.append((vec[i],i))
    values = np.sort(np.array(values,dtype),order='val') 
    new_points = []
    new_ind = []
    for tup in values:
        new_points.append(points[tup[1]])
        new_ind.append(ind[tup[1]])
    return np.array(new_points),np.array(new_ind)




def M(points,tri,Nint):
    
    indptr,ind = tri.vertex_neighbor_vertices
    W = np.zeros((Nint,Nint)) # 余切权值矩阵
    A = np.zeros((Nint,1)) # 每个顶点i对应的贡献面积数组A[i]
    
    for i in range(Nint):
        tot = 0
        nhb_ind = ind[indptr[i]:indptr[i+1]] # 顶点i的邻点索引列表
        nhb = points[nhb_ind] # 邻点坐标
        nhb,nhb_ind = tri_ang(nhb,nhb_ind,points[i]) # 对i的所有邻点按极角排序
        
        for j in range(len(nhb_ind)):
            vec = nhb[j]-points[i] # 顶点i到排序后第j个邻点的向量
            vec_av = nhb[j-1]-points[i] # 顶点i到排序后第j-1个邻点的向量
            if j+1 >= len(nhb_ind):
                k = 0
            else:
                k = j+1
            vec_ap = nhb[k]-points[i] # 顶点i到排序后第j+1个邻点的向量
            
            A[i] += 0.5/3*np.linalg.norm(np.cross(vec,vec_av))
            
            if nhb_ind[j] < Nint:
                # 通过点积和叉积计算余切值：A.B/||AxB||
                cotan_alpha = np.dot(vec_av,vec_av-vec)/np.linalg.norm(np.cross(vec_av,vec_av-vec))
                cotan_beta = np.dot(vec_ap,vec_ap-vec)/np.linalg.norm(np.cross(vec_ap,vec_ap-vec))
                # 计算权值矩阵元素Wij
                W[i,nhb_ind[j]] = -0.5*(cotan_alpha+cotan_beta)
                
                tot += cotan_alpha+cotan_beta
    
        W[i,i] = -0.5*tot # 计算矩阵对角线元素
    
    return (1/A)*W


def rect(x,y,L,H,x0=0,y0=0):
    if 0<x-x0<L and 0<y-y0<H:
        return -1
    else:
        return 1

def rect_rim(L,H,a,x0=0,y0=0):
    tab1 = np.arange(x0,L+x0,a)[:,np.newaxis]
    h = np.hstack((tab1,H*np.ones((len(tab1),1))+y0))
    b = np.hstack((tab1,np.zeros((len(tab1),1))+y0))
    tab2 = np.arange(y0+a,H+y0,a)[:,np.newaxis]
    g = np.hstack((np.zeros((len(tab2),1))+x0,tab2))
    d = np.hstack((L*np.ones((len(tab2),1))+x0,tab2))
    hp = np.array([[L+x0,H+y0]])
    bp = np.array([[L+x0,0]])
    return np.vstack((h,b,g,d,hp,bp))
    


# 1*1方形区域测试示例
L = 1
H = 1

dl = 0.05
sol = in_curve(rect,[L,H],(100,100),dl)
sol_tri,Nint,Nbord = triang(sol,dl,rect,[L,H],rect_rim(L,H,dl))

# 构建三角剖分
tri = scipy.spatial.Delaunay(sol_tri)

M = M(sol_tri,tri,Nint)

valp,vecp = np.linalg.eig(M) # 求解特征值与特征向量
vecp = np.real(vecp)

# 与解析解对比
T = 1000
U = np.arange(0,T,1)
NUsim = np.array([len(count_vp(valp,u)) for u in U])
NU = np.array([len(rect_drum(L,H,u)) for u in U])
plt.plot(U,NUsim,label='simulation')
plt.plot(U,NU,label='exacts')
plt.legend()
plt.show()


# 某阶特征向量的3D可视化
vecp_tot = np.vstack((vecp,np.zeros((Nbord,Nint))))
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(sol_tri[:,0],sol_tri[:,1],vecp_tot[:,0],triangles=tri.simplices)
plt.show()

内容的提问来源于stack exchange，提问作者bosonic