# Applying Neumann boundaries to Crank-Nicolson solution in python

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
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 = B+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())

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 = 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.

• is there a typo in the first set of BC (missed $_x$) for the $u(l,t)$? – Anton Menshov Mar 24 '18 at 16:25
• @AntonMenshov no, it's $\partial u/\partial x = 0$ at $x=l$ – Googlebot Mar 24 '18 at 16:29
• You have both $u(l,t)=0$ and $u_x(l,t)=0$ at the same time? That sounds strange. – Anton Menshov Mar 24 '18 at 16:32
• @AntonMenshov the former is the boundary condition, which has been used in the code. I want to replace it with the latter condition to modify the code. – Googlebot Mar 24 '18 at 16:34
• Oh got it. At this moment, you have Dirichlet $u(l,t)=0$ and you want to change it to Neumann $u_x(l,t)=0$. Becuase the combination is also possible - that's Robin BC. – Anton Menshov Mar 24 '18 at 16:35

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}.$$.