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:
    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:

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.

  • $\begingroup$ is there a typo in the first set of BC (missed $_x$) for the $u(l,t)$? $\endgroup$
    – Anton Menshov
    Commented Mar 24, 2018 at 16:25
  • $\begingroup$ @AntonMenshov no, it's $\partial u/\partial x = 0$ at $x=l$ $\endgroup$
    – Googlebot
    Commented Mar 24, 2018 at 16:29
  • $\begingroup$ You have both $u(l,t)=0$ and $u_x(l,t)=0$ at the same time? That sounds strange. $\endgroup$
    – Anton Menshov
    Commented Mar 24, 2018 at 16:32
  • $\begingroup$ @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. $\endgroup$
    – Googlebot
    Commented Mar 24, 2018 at 16:34
  • $\begingroup$ 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. $\endgroup$
    – Anton Menshov
    Commented Mar 24, 2018 at 16:35

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


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.