Python网格二重积分（带φ(y)dxdy测度）稳定性及PDE数值发散问题问询

给定网格$((x_i,y_j)){i,j}$，函数$f$在网格点处取值已知，需计算$f$关于测度$\phi(y)dxdy$的二重积分（$\phi$为归一化函数，$dx=x{i+1}-x_i$，$dy=y_{j+1}-y_j$）。

已实现矩形法积分函数如下：

def integral(f):
    S=0
    for i in range(len(x)):
        for j in range(len(y)):
            S+=f[i,j]*phi[j]*dx*dy
    return S

该积分函数是PDE求解的核心模块（需调用百次以上），但稳定性不佳。涉及的PDE为带边界条件的输运方程，边界由依赖解的两个二重积分乘积给出。

测试示例：$x、y$非负，$f(x,y)=2e^{{-(x+y)}$，$\phi(y)=e}{-y}$。

PDE数值格式与初始/边界条件 #

数值格式为：
$$u^{{n+1}_{i,j}=\frac{dt}{dx}u}n_{i-1,j}+\left(1-\frac{dt}{dx}\right)u^{n_{i,j}-\gamma(x_i,y_j)F(t_n)u}n_{i,j}dt$$
其中：

• 边界条件：$u^n_{0,j}=F(t_n)$
• 初始条件：$u^0_{i,j}=f(x_i,y_j)$
• $F(t_n)=\int_{\mathbb{R}^{+\times\mathbb{R}}+}g(x,y)u(t_n,x,y)dxdy$

当前问题 #

单步计算结果正常，但迭代求解PDE时数值随时间急剧增大，寻求解决方案（注：solve_pde函数中仅随时间更新网格上的$f$值）。

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