# High frequency noise at solving diffusion equation

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[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)


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.

• Have you tried decreasing your timestep? It looks like you're probably violating the relevant CFL condition. – Doug Lipinski Apr 23 '14 at 15:46
• Yes, this actually solves the problem. Thank's a lot. – sonium Apr 23 '14 at 15:54
• Glad to help, I'll post this up as an answer. – Doug Lipinski Apr 23 '14 at 15:54

In the case of the heat equation (Fick's 2nd law is a heat equation), your CFL number $C$ must satisfy
$$C= k\dfrac{\Delta t}{(\Delta x)^2}<A$$
where $k$ is the diffusion coefficient and $A$ is some constant that is typically in the range (0,1] and depends on the specific equation and method used. This places a somewhat harsh constraint on your timestep as you increase the number of grid points (decreasing $\Delta x$) and is one of the major reasons implicit methods are used.