High frequency noise at solving differential equation

Question

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.

Roth

Answer 1

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

The easiest way to explain this is for the 1-D case: assume your concentration is a function of spatial coordinate x and timestep i. If you do all the math (write down your equations, substitute the partial derivatives with finite differences, 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)]

so the concentration of a point on the next step depends on its previous value and the ones of its two neighbors. It is not too hard to see that when k = 0.5, every point gets replaced by the average of its two neighbors, so a concentration profile of [...,0,1,0,1,0,...] will become [...,1,0,1,0,1,...] on the next step. If k > 0.5, such a profile will blow up exponentially. You calculate your second order derivative with a longer convolution (I effectively use [1,-2,1]), but I guess that does not change anything for the instability problem.

I don't know about normal diffusion, but based on experience with thermal diffusion, I would guess that k scales with dt * diffusion_coeff / dx^2. You thus have to chose your timestep small enough so that your simulation does not become instable. To make the simulation stable, but still as fast as possible, chose your parameters so that k is a bit smaller than 0.5. Something similar can be derived for 2-D and 3-D cases. The easiest way to achieve this is to increase dx, since your total calculation time will scale with 1/dx^3 for a linear problem, 1/dx^4 for 2-D problems, and even 1/dx^5 for 3-D problems.

There are better methods to solve diffusion equations, I believe that Crank Nicolson is at least standard for solving heat-equations (which is also a diffusion problem). The 'problem' is that this is an implicit method, which means that you have to solve a set of equations to calculate your 'concentration' at the next timestep, which is a bit of a pain to implement. But this method is guaranteed to be numerical stable, even for big timesteps.

Landgravine

Answer 2

Using scipy.solve_ivp to solve this initial value problem will use a more sophisticated algorithm, thereby avoiding the cumulative errors that lead to a nonsensical solution. I have also added necessary boundary conditions that imposes a flux of zero at both extremities.

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from scipy.interpolate import splrep, splev
gridpoints = 128

def profile(x):
    peak_x = 2.
    straggle = .1576
    dose = 1 
    return dose/(np.sqrt(2*np.pi)*straggle)*np.exp(-(x-peak_x)**2/2/straggle**2)

def numpy_dy_dx(y, x):
    return np.gradient(y, x[1] - x[0])

def spline_dy_dx(y, x):
    tck = splrep(x, y, k=2)
    return splev(x, tck, der=1)

def dydt(t, y, x, diffusion_coefficient):
    dy_dx = spline_dy_dx(y, x)
    # Boundary conditions dy/dx at both extremities = 0.0
    dy_dx[0] = 0.0
    dy_dx[-1] = 0.0
    dy2_dx2 = spline_dy_dx(dy_dx, x)
    dy_dot = diffusion_coefficient * dy2_dx2
    return dy_dot
    
x = np.linspace(0, 4, gridpoints)
nx = profile(x)

plt.figure(figsize=(12,8))

steps = 10
diffusion_coefficient = 0.002
t0 = 0.0 
tf = 500.0
t_eval = np.linspace(t0, tf, steps)
t_span = (t0, tf)

solution = solve_ivp(dydt, t_span, nx, t_eval=t_eval, args=(x, diffusion_coefficient), method='Radau')
print(f'Number of dydt function evaluations: {solution.nfev}')

for i in range(len(solution.t)):
    plt.plot(x, solution.y[:, i], label=f't = {solution.t[i]:0.0f}')
    
plt.grid()
plt.legend()

# Verify that energy is conserved
from scipy.integrate import simpson
v1 = simpson(nx, dx=x[1] - x[0])
for i in range(len(solution.t)):
    v2 = simpson(solution.y[:, i], dx=x[1] - x[0])
    print(solution.t[i], v2)

The following is the resulting plot. The solver used 149 function evaluations:

Note: derivatives can either be calculated with numpy_dy_dx which makes use of numpy.gradient function or with spline_dy_dx which makes use of scipy.interpolate.splrep and scipy.interpolate.splev.

Stuff