# Numerical artefacts in solution of spherical heat equation using FDM

I was attempting to solve the diffusion equation for a solid sphere using a naive FDM scheme. The governing PDE for the scalar concentration field $$u(r,t)$$ is $$u_t = r^{-2}(r^2 \alpha u_r)_r, \quad r \in [0,1],$$ and is accompanied by the boundary conditions $$u_r(0,t) = 0, \quad u(1,t) = 0$$ The initial condition is $$u(r,1) = 1$$. For simplicity $$\alpha$$ can be assumed to be constant.

I implemented the scheme given in the book of Langtangen and Linge (2017), Finite Difference Computing with PDEs, Springer. An exemplary Python program is given below:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def sphere_fdm(t,y,r,alpha,dr):
"""Naive FDM discretization for spherically-symmetric diffusion

References:
[1] Langtangen, H. P., & Linge, S. (2017). Finite difference computing
with PDEs: a modern software approach. Springer Nature. (pg. 252)
"""
n = y.shape[0]
dydt = np.empty_like(y)

# Symmetry condition
dydt[0] = 6*alpha*(y[1] - y[0])/dr**2

for i in range(1,n-1):

rl = 0.5*(r[i] + r[i-1])
rr = 0.5*(r[i+1] + r[i])

gl = rl**2 * alpha * (y[i] - y[i-1])
gr = rr**2 * alpha * (y[i+1] - y[i])

dydt[i] = (1./r[i])**2 * (gr - gl)/dr**2

# Dirichlet boundary
dydt[-1] = 0

return dydt

def sphere_sol(t,r,alpha,nterms=40):
"""Series solution of heat equation with Dirichlet boundaries

References:
[1] Crank, J. (1979). The Mathematics of Diffusion.
Oxford University Press. (pg. 91)
"""
a_ = 1
C = np.zeros_like(r)
csum = 0
for n in range(1,nterms+1):
expterm = np.exp(-alpha*n**2*np.pi**2*t/a_**2)
C[1:] += (-1)**n/n * np.sin(n*np.pi*r[1:]/a_) * expterm
csum += (-1)**n * expterm

C[1:] *= 2*a_/(np.pi*r[1:])
C[0] = 2*csum

C += 1
return 1 - C

def main():

N = 41
r, dr = np.linspace(0,1,N,retstep=True)

t = 0.1
alpha = 0.1
method="RK45"

y = np.ones(N)
y[-1] = 0

rhs = lambda t, y: sphere_fdm(t,y,r,alpha,dr)
sol = solve_ivp(rhs,[0,t],y,method=method)

fig, ax = plt.subplots()
iax = ax.inset_axes([0.1,0.2,0.5,0.5])

ax.plot(r,sol.y)
ax.plot(r,sphere_sol(t,r,alpha),'o')

ax.set_xlabel("r")
ax.set_ylabel("Y(r,t)")

iax.plot(r,sol.y)
iax.set_xlim(0,3*dr)
iax.set_ylim(0.98,1.02)
iax.set_xlabel("r")

plt.show()

if __name__ == '__main__':
main()


For the solution in time, the method of lines (MOL) is used with the solvers from scipy.integrate.solve_ivp.

While the solution agrees fairly well with the analytical solution (the red dots), when using the "RK45" integrator, the solution develops an unexpected kink at the center:

The kink is not observable when I use the "BDF" method:

Is there any explanation for this behavior (and assuming it's not an error in my discretized equations or Python code)?

• It's hard to tell what's happening temporally. Does the kink develop over time or does it occur at early time then go away? May 30, 2023 at 21:29
• The kink appears gradually in the penetration period (short diffusion times) and then slowly disappears in the long-term regime. I suspect it may be an artefact from the jump in the initial condition, exacerbated by the singularity near the origin. Langtangen mentions that Rannacher time-stepping can eliminate such problems, similar to the "BDF" method in my case. Other threads on SE suggest implicit methods being more suited to parabolic PDEs. May 30, 2023 at 22:05

You may simply change


def sphere_fdm(t,y,r,alpha,dr):

# Symmetry condition
dydt[0] = 6*alpha*(y[1] - y[0])/dr**2

# Loop
for i in range(1,n-1):
...

# Dirichlet boundary
dydt[-1] = 0

def main():

y = np.ones(N)
y[-1] = 0



to


def sphere_fdm(t,y,r,alpha,dr):

# Loop
for i in range(1,n-1):
...

# Symmetry condition
dydt[0] = dydt[1]

# Dirichlet boundary
dydt[-1] = 0

def main():

y = np.ones(N)
y[0]  = y[1]
y[-1] = 0



RK45

BDF