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


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,
                    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

output with noise

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.

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

Solving partial differential equations with explicit timestepping methods relies on meeting a certain CFL condition for stability. Since you are using a forward Euler timestepping scheme, you must ensure that your timestep is small enough compared to the grid spacing.

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.

Note that this condition is necessary for stability, meaning that the numerical solution will eventually blow up if this condition is not satisfied. Concerns about accuracy are another topic, but if the solution has blown up you're guaranteed to be inaccurate.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.