Application of Boundary Conditions in finite difference solution for the heat equation and Crank-Nicholson

Question

The code below solves the 1D heat equation that represents a rod whose ends are kept at zero temparature with initial condition 10*np.sin(np.pi*x).

How are the Dirichlet boundary conditions (zero temp at both ends) worked into this calculation? I was told the upper, lower rows of matrix A contain two non-zero elements, and the missing third element is the Dirichlet condition. But I do not understand by which mechanism this condition affects the calculation. With missing elements in A, how can u_{0} or u_{n} be zero?

The finite difference method below uses Crank-Nicholson.

import numpy as np
import scipy.linalg

# Number of internal points
N = 200

# Calculate Spatial Step-Size
h = 1/(N+1.0)
k = h/2

x = np.linspace(0,1,N+2)
x = x[1:-1] # get rid of the '0' and '1' at each end

# Initial Conditions
u = np.transpose(np.mat(10*np.sin(np.pi*x)))

# second derivative matrix
I2 = -2*np.eye(N)
E = np.diag(np.ones((N-1)), k=1)
D2 = (I2 + E + E.T)/(h**2)

I = np.eye(N)

TFinal = 1
NumOfTimeSteps = int(TFinal/k)

for i in range(NumOfTimeSteps):
    # Solve the System: (I - k/2*D2) u_new = (I + k/2*D2)*u_old
    A = (I - k/2*D2)
    b = np.dot((I + k/2*D2), u)
    u = scipy.linalg.solve(A, b)

Answer 1

Let's have a look at a simple example. We assume N = 3, i.e. three inner points, but we will first also include the boundary points in the matrix D2 describing the approximate second derivatives:

      1  /  1 -2  1  0  0 \
D2 = --- |  0  1 -2  1  0 |
     h^2 \  0  0  1 -2  1 /

The first line means the approximate second derivative at x_1 is 1/h^2 * (u_0 - 2*u_1 + u_2). We know that u_0 = 0 though, due to the homgeneous Dirichlet boundary conditions, so we can simply omit it from the equation, and e get the same result for the matrix

      1  /  0 -2  1  0  0 \
D2 = --- |  0  1 -2  1  0 |
     h^2 \  0  0  1 -2  0 /

Since u_0 and u_{n+1} are not real unknowns -- they are known to be zero -- we can completely drop them from the matrix, and we get

      1  /  2  1  0 \
D2 = --- |  1 -2  1 |
     h^2 \  0  1 -2 /

The missing entries in the matrix really correspond to the fact that the boundary conditions are zero.