I'm trying to simulate a simple diffusion based on 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 - x # 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)
for the first few time steps everything looks fine, but then I start to get high frequency noise, do to build-up from numerical errors which are amplified through the second derivative. Since it seems to be hard to increase the float precision I'm hoping that there is something else that I can do to suppress this? I already increased the number of points that are being used to construct the 2nd derivative.