6
$\begingroup$

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: Spherical diffusion with RK45

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

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

$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – whpowell96
    May 30, 2023 at 21:29
  • $\begingroup$ 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. $\endgroup$
    – IPribec
    May 30, 2023 at 22:05

1 Answer 1

2
$\begingroup$

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

RK45

BDF

BDF

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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