求Python中与Matlab del2等效的离散拉普拉斯计算方法

Matlab的del2函数用于计算离散拉普拉斯算子，核心逻辑是基于邻域点加权平均与中心点的差值，不同维度的计算公式如下：

• 1D数组：del2(u) = (u[i-1] + u[i+1] - 2*u[i]) / (2*h²)（h为步长，默认h=1）
• 2D数组：del2(u) = (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1] - 4*u[i,j]) / (4*h²)（默认h=1）

下面分两种场景给出实现方案：

一、用NumPy实现通用的del2等效函数 #

如果处理常规数值数组，直接用NumPy就能实现和del2完全一致的功能：

import numpy as np

def del2(u, h=1.0):
    # 处理1D数组
    if u.ndim == 1:
        n = len(u)
        laplacian = np.zeros_like(u)
        # 内部点计算
        laplacian[1:-1] = (u[:-2] + u[2:] - 2*u[1:-1]) / (2 * h**2)
        # 边界点（对齐Matlab默认的一阶差分近似）
        laplacian[0] = (u[1] - u[0]) / h**2
        laplacian[-1] = (u[-2] - u[-1]) / h**2
        return laplacian
    # 处理2D数组
    elif u.ndim == 2:
        rows, cols = u.shape
        laplacian = np.zeros_like(u)
        # 内部点计算
        laplacian[1:-1, 1:-1] = (u[:-2,1:-1] + u[2:,1:-1] + u[1:-1,:-2] + u[1:-1,2:] - 4*u[1:-1,1:-1]) / (4 * h**2)
        # 上下边界点
        laplacian[0, 1:-1] = (u[1,1:-1] - u[0,1:-1]) / h**2
        laplacian[-1, 1:-1] = (u[-2,1:-1] - u[-1,1:-1]) / h**2
        # 左右边界点
        laplacian[1:-1, 0] = (u[1:-1,1] - u[1:-1,0]) / h**2
        laplacian[1:-1, -1] = (u[1:-1,-2] - u[1:-1,-1]) / h**2
        # 四个角落点
        laplacian[0,0] = (u[0,1] + u[1,0] - 2*u[0,0]) / (2 * h**2)
        laplacian[0,-1] = (u[0,-2] + u[1,-1] - 2*u[0,-1]) / (2 * h**2)
        laplacian[-1,0] = (u[-1,1] + u[-2,0] - 2*u[-1,0]) / (2 * h**2)
        laplacian[-1,-1] = (u[-1,-2] + u[-2,-1] - 2*u[-1,-1]) / (2 * h**2)
        return laplacian
    else:
        raise ValueError("仅支持1D和2D数组")

使用示例：

# 1D测试
u1d = np.array([0, 1, 0])
print(del2(u1d))  # 输出: [-0.5 -1.  -0.5]，与Matlab del2(u1d)结果一致

# 2D测试
u2d = np.array([[0,0,0],[0,1,0],[0,0,0]])
print(del2(u2d))  # 中心值为-1，边界角落为-0.5，其余边界为-0.5，与Matlab结果一致

二、在DeepXDE框架中的实现（结合你的代码示例） #

你的代码基于DeepXDE求解耦合PDE，已用到dde.grad.hessian计算二阶导数。在该框架下实现del2等效逻辑，需结合自动微分的符号计算特性：

1. 1D场景（对应你的代码） #

1D中拉普拉斯算子即二阶导数d²y/dx²，Matlab的del2(y)等于该二阶导数的1/2（h=1时）。只需对已有的二阶导数做系数调整：

import deepxde as dde

# 假设已定义参数：D, k, a, epsilon, mu_1, mu_2, b

def pde_1D(x, y):
    V, W = y[:, 0:1], y[:, 1:2]
    dv_dt = dde.grad.jacobian(y, x, i=0, j=1)
    # 计算V的二阶导数（拉普拉斯算子）
    laplacian_V = dde.grad.hessian(y, x, component=0, i=0, j=0)
    # 等效于Matlab del2(V)（h=1时）
    del2_V = laplacian_V / 2
    dw_dt = dde.grad.jacobian(y, x, i=1, j=1)
    
    ## 耦合PDE+ODE方程
    # 若原方程需用del2结果，替换dv_dxx为del2_V即可
    eq_a = dv_dt - D * del2_V + k*V*(V-a)*(V-1) + W*V
    eq_b = dw_dt - (epsilon + (mu_1*W)/(mu_2+V))*(-W -k*V*(V-b-1))
    return [eq_a, eq_b]

2. 多维场景 #

2D/3D中拉普拉斯算子是各维度二阶导数的和，del2等效值为拉普拉斯算子除以对应维度系数（2D除以4，3D除以6，h=1时）。例如2D场景：

# 2D下计算del2等效值
laplacian_V = dde.grad.hessian(y, x, component=0, i=0, j=0) + dde.grad.hessian(y, x, component=0, i=1, j=1)
del2_V = laplacian_V / 4

注意事项 #

• 当步长h≠1时，需在系数中加入h²分母，例如1D：del2_V = laplacian_V / (2*h**2)，2D：del2_V = laplacian_V / (4*h**2)。
• NumPy版本严格对齐Matlabdel2的边界处理逻辑，若无需边界近似，可仅计算内部点。

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