High-frequency noise when solving a differential equation

Question

I am trying to simulate a simple diffusion based on Fick 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) 

during the first few steps of the time everything looks fine, but then I start to receive high-frequency noise, to do it in order to increase from numerical errors, which are amplified through the second derivative. Since it seems hard to increase the accuracy of the float, I hope I can do something else to suppress this? I have already increased the number of points that are used to construct the 2nd derivative.

Answer 1

I don't have time to study your solution in detail, but it looks like you are solving a partial differential equation with a forward Euler scheme . This is pretty easy to implement, as you show, but it can become numerically unstable if your timestep is too small. Your only solution is to reduce the time interval or increase the spatial resolution.

The easiest way to explain this for the 1-D case: suppose your concentration is a function of the spatial coordinate x and timestep i . If you do all the math (write down your equations, replace partial derivatives with finite differences , it should be pretty easy), you will probably get something like this:

 C(x, i+1) = [1 - 2 * k] * C(x, i) + k * [C(x - 1, i) + C(x + 1, i)] 

therefore, the concentration of a point in the next step depends on its previous value and the values ​​of its two neighbors. It is easy to see that at k = 0.5 each point is replaced by the average number of its two neighbors, therefore, in the next step, the concentration profile [...,0,1,0,1,0,...] becomes [...,1,0,1,0,1,...] . If k > 0.5 , such a profile will explode exponentially. You compute a second-order derivative with a longer convolution (I use [1, -2,1] effectively), but I think this does not change anything for the instability problem.

I don’t know about normal diffusion, but based on experience with thermal diffusion, I would suggest that k scales with dt * diffusion_coeff / dx^2 . Therefore, you must choose your own time period small enough so that your simulation does not become unstable. To make the simulation stable, but still as fast as possible, you have chosen your options so that k slightly less than 0.5 . Something similar can be obtained for 2-D and 3-D cases. The easiest way to achieve this is to increase dx , since your total calculation time will be scaled with 1/dx^3 for a linear problem, 1/dx^4 for two-dimensional problems, and even 1/dx^5 for three-dimensional problems.

There are more efficient methods for solving diffusion equations, I believe that Crank Nicolson is at least standard for solving heat equations (which is also diffusion). The “problem” is that it is an implicit method, which means that you need to solve a set of equations to calculate your “concentration” for the next time period, which is a little painful to implement. But this method is guaranteed to be numerically stable even for large time intervals.