求解微分方程时的高频噪声

Question

我正在尝试模拟一个简单的扩散过程，这个过程是基于菲克第二定律。

from pylab import *
import numpy as np
gridpoints = 128

def profile(x):
    range = 2.
    straggle = .1576
    dose = 1 
    return dose/(sqrt(2*pi)*straggle)*exp(-(x-range)**2/2/straggle**2)

x = linspace(0,4,gridpoints)
nx = profile(x)
dx = x[1] - x[0] # use np.diff(x) if x is not uniform
dxdx = dx**2

figure(figsize=(12,8))

plot(x,nx)
timestep = 0.5
steps = 21
diffusion_coefficient = 0.002
for i in range(steps):
    coefficients = [-1.785714e-3, 2.539683e-2, -0.2e0, 1.6e0,
                    -2.847222e0,
                    1.6e0, -0.2e0, 2.539683e-2, -1.785714e-3]
    ccf = (np.convolve(nx, coefficients) / dxdx)[4:-4] # second order derivative
    nx = timestep*diffusion_coefficient*ccf + nx
    plot(x,nx)

在前几个时间步骤中，一切看起来都很好，但之后我开始出现高频噪声，这是由于数值误差的累积，这些误差通过二阶导数被放大。因为提高浮点数的精度似乎很困难，所以我希望能找到其他方法来抑制这种噪声。我已经增加了用于构建二阶导数的点的数量。

Answer 1

我没有时间详细研究你的解决方案，但看起来你是在用前向欧拉法来解偏微分方程。这种方法实现起来比较简单，正如你所展示的，但如果你的时间步长太小，可能会导致数值不稳定。你唯一的解决办法就是减小时间步长或者增加空间分辨率。

最简单的解释是针对一维情况：假设你的浓度是空间坐标x和时间步i的函数。如果你把所有的数学运算都写下来（把方程写出来，用有限差分法替代偏导数，这应该很简单），你可能会得到类似这样的结果：

C(x, i+1) = [1 - 2 * k] * C(x, i) + k * [C(x - 1, i) + C(x + 1, i)]

所以，下一步某个点的浓度依赖于它之前的值和它两个邻居的值。很容易看出，当k = 0.5时，每个点的值会被它两个邻居的平均值替代，因此浓度分布[...,0,1,0,1,0,...]在下一步会变成[...,1,0,1,0,1,...]。如果k > 0.5，这样的分布会指数级增长。你用更长的卷积来计算二阶导数（我实际上使用的是[1,-2,1]），但我想这对不稳定性问题没有改变什么。

我对正常扩散不太了解，但根据我在热扩散方面的经验，我猜k与dt * diffusion_coeff / dx^2成比例。因此，你必须选择一个足够小的时间步长，以确保你的模拟不会变得不稳定。为了让模拟稳定，同时尽可能快，选择你的参数使得k稍微小于0.5。对于二维和三维情况也可以推导出类似的结论。实现这一点最简单的方法是增加dx，因为对于线性问题，你的总计算时间将与1/dx^3成比例，对于二维问题是1/dx^4，而对于三维问题则是1/dx^5。

我相信还有更好的方法来解决扩散方程，我认为Crank Nicolson方法至少是解决热方程（这也是一个扩散问题）的标准方法。问题在于这是一种隐式方法，这意味着你需要解一组方程来计算下一时间步的“浓度”，这实现起来有点麻烦。但这种方法可以保证数值稳定，即使在较大的时间步长下也是如此。