有限差分法与Crank-Nicholson法中边界条件的应用

Question

下面的代码解决了一个一维热方程，这个方程描述的是一根两端温度保持在零度的杆，初始条件是10*np.sin(np.pi*x)。

那么，怎么把边界条件（两端温度为零）融入到这个计算中呢？有人告诉我，矩阵A的上下两行包含两个非零元素，而缺失的第三个元素就是这个边界条件。但我不明白这个条件是通过什么机制影响计算的。既然A中缺少元素，u_{0}或u_{n}怎么能是零呢？

下面的有限差分方法使用了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

有人告诉我，矩阵A的上面和下面的行包含两个非零元素，而缺失的第三个元素（也就是零）是Dirichlet条件。

我不太确定你听到的是什么，但我可以告诉你我对这个问题的理解。

Dirichlet边界条件是用来固定某个值的，比如在这个例子中是温度。Neumann边界条件则是指定某个点的流量或一阶导数。如果你在某个表面有对流边界条件，你就需要用到这个。

关于如何处理Dirichlet边界条件，首先你要在不考虑边界条件的情况下来构建系统矩阵。如果你在某个节点上有一个固定的温度，你可以这样处理：

1. 把右侧向量中对应节点的行设置为边界温度。然后把左侧矩阵中该行的所有列都设为零，并把对应节点的对角元素设为1。

Answer 2

让我们来看一个简单的例子。我们假设 N = 3，也就是有三个内部点，但我们首先还会在描述近似二阶导数的矩阵 D2 中包含边界点：

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

第一行的意思是，在 x_1 处的近似二阶导数是 1/h^2 * (u_0 - 2*u_1 + u_2)。不过，由于我们有均匀的 Dirichlet 边界条件，我们知道 u_0 = 0，所以可以把它从这个公式中省略掉，这样我们得到的矩阵结果是一样的。

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

因为 u_0 和 u_{n+1} 不是我们需要求解的未知数——它们的值已经知道是零——所以我们可以完全把它们从矩阵中去掉，这样我们就得到了：

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

矩阵中缺失的部分实际上对应于边界条件为零的事实。