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)

Garmon

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.

Ushas

Then would it be correct to say the code removes two values from beginning and the end of x, with x = x[1:-1] for this reason? So that we get a new "reduced" u that can multiply the "reduced" A? – Garmon

@user423805 - yes. And if the imposed boundary conditions were of a form where you couldn't so easily construct a reduced A, then you'd have to consider the full A and x. In this (quite common) case, the values of x corresponding not to real data but of ways of applying the boundary condition are commonly called ghost cells or guard cells or sentinal cells (or nodes, or zones, or what have you). – Laverty

Answer 2

I was told the upper, lower rows of matrix A contain two non-zero elements, and the missing third element (that is zero) is the Dirichlet condition.

I'm not sure what you've been told, but here's what I know about this problem.

Dirichlet boundary conditions fix the value of the potential (temperature in this case). A Neumann boundary condition will specify flux or first derivative at a point. You'll need this if you have convection boundary conditions at a surface.

As for treating Dirichlet boundary conditions, you formulate the system matrix without considering boundary conditions first. If you have a fixed temperature at a given node, you can handle it this way:

1. Set the row in the right hand side vector for the given node equal to the boundary temperature. Zero out all columns of that row in the left hand side matrix and set the diagonal element that corresponds to that node equal to one.

Riegel