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

Question

我试图模拟一个基于Fick's 2nd law的简单扩散。在

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

我没有时间详细研究你的解，但你似乎在用forward Euler scheme来解偏微分方程。如您所示，这是非常容易实现的，但是如果时间步长太小，这可能会在数值上变得不稳定。您唯一的解决方案是减少时间步长或提高空间分辨率。在

解释这一点最简单的方法是一维情况：假设你的注意力是空间坐标x和时间步长i的函数。如果你做了所有的数学运算（写下你的方程，用finite differences代替偏导数，应该很容易），你可能会得到这样的结果：

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,...]的浓度分布在下一步将变成{}。如果k > 0.5，这样的配置文件将以指数方式爆炸。你用一个较长的卷积计算二阶导数（我有效地使用了[1，-2,1]），但我想这对不稳定性问题没有任何改变。在

我不知道正常扩散，但根据热扩散的经验，我猜k与{}成比例。因此，您必须选择足够小的时间步长，以便模拟不会变得不稳定。为了使模拟稳定，但仍然尽可能快，请选择参数，使k比0.5小一点。在二维和三维情况下也可以得出类似的结论。实现这一点最简单的方法是增加dx，因为对于线性问题，您的总计算时间将按1/dx^3来调整，1/dx^4对于二维问题，甚至对于三维问题，1/dx^5。在

有更好的方法来求解扩散方程，我相信Crank Nicolson至少是求解热方程的标准（这也是一个扩散问题）。问题是这是一个implicit方法，这意味着你必须求解一组方程来计算下一个时间步的“浓度”，这有点难以实现。但这种方法保证了数值稳定性，即使在大时间步长下也是如此。在