Applying Neumann boundaries to Crank-Nicolson solution in python

Question

Consider the heat equation $$u_t = \kappa u_{xx}$$ with boundary conditions of $$u(x,0)=0\\ u(0,t)=100\\ u(l,t)=0$$ Numerical analysis by pyton can be done with

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags

def Crank_Nicolson(dy,ny,dt,nt,k,T,ntout):
    Tout = []
    T0 = T[0]
    T1 = T[-1]
    s = k*dt/dy**2
    A = diags([-0.5*s, 1+s, -0.5*s], [-1, 0, 1], shape=(ny-2, ny-2)).toarray() 
    B1 = diags([0.5*s, 1-s, 0.5*s],[-1, 0, 1], shape=(ny-2, ny-2)).toarray()
    for n in range(1,nt):
        Tn = T
        B = np.dot(Tn[1:-1],B1) 
        B[0] = B[0]+0.5*s*(T0+T0)
        B[-1] = B[-1]+0.5*s*(T1+T1)
        T[1:-1] = np.linalg.solve(A,B)
        if n % int(nt/float(ntout)) == 0 or n==nt-1:
            Tout.append(T.copy())
    return Tout,s

dt = 0.01
dy = 0.001
k = 10**(-4)
y_max = .1
y = np.arange(0,y_max+dy,dy) 
ny = len(y)
nt = 1000

T = np.zeros((ny,))
T[0] = 100
Tout,s = Crank_Nicolson(dy,ny,dt,nt,k,T,10)

for T in Tout:
    plt.plot(y,T)    
plt.show()

How a Neuman boundary such as $$u_x(l,t)=0$$ can be implemented in this code?

This a case of keeping one end of a bar at temperature 100 °C while the other end is insulated.

Answer 1

Boundary conditions are often an annoyance, and can frequently take up a surprisingly large percentage of a numerical code. How to implement them depends on your choice of numerical method. Finite difference schemes often find Dirichlet conditions more natural than Neumann ones, whereas the opposite is often true for finite element and finite methods applied to diffusive problems.

One standard technique, which will work well in your case is to introduce an additional conceptually virtual degree of freedom or "ghost" point lying at $x=l+\Delta x$, then apply the model equation at $x = l$ in the form $$ \frac{du(l)}{dt}=\kappa(l)\frac{u(l+\Delta x)-2u(l)+u(l-\Delta x)}{{\Delta x}^2}$$ as well as boundary condition at $x = l$, in the form $$\frac{u(l+\Delta x,t)-u(l-\Delta x,t)}{2\Delta x} = 0.$$

Note that there is no explicit equation applied at $x=l+\Delta x$, so there are still as many equations as unknowns, so the system remains well posed. In fact, most practical implementations rearrange the boundary condition into the form $$u(l+\Delta x)=u(l-\Delta x),$$ and just substitute back into the previous equation to get $$ \frac{du(l)}{dt}=\kappa(l)\frac{2u(l-\Delta x)-2u(l)}{{\Delta x}^2}.$$.