我正在尝试基于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)

在最初的几个时间步骤中，一切看起来都很好，但随后我开始得到高频噪声，这是从通过二阶导数放大的数值误差中积累起来的。由于似乎很难提高浮点精度，我希望还有其他方法可以抑制这种情况？我已经增加了用于构造二阶导数的点数。