Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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)

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.

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

1 Answer 1

up vote 8 down vote accepted

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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