6

下面的代码求解一维热方程，该方程表示一根杆，其末端保持在零温度，初始条件为 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)
4

2 回答 2

5

让我们看一个简单的例子。我们假设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)。我们知道u_0 = 0，由于齐次狄利克雷边界条件，所以我们可以简单地从方程中省略它，并且 e 对于矩阵得到相同的结果

      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 /

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

于 2011-01-30T13:39:32.247 回答
0

有人告诉我，矩阵 A 的上下两行包含两个非零元素，而缺少的第三个元素（即零）是狄利克雷条件。

我不确定你被告知了什么，但这是我对这个问题的了解。

Dirichlet 边界条件确定了势的值（在这种情况下为温度）。Neumann 边界条件将指定一个点的通量或一阶导数。如果您在表面上有对流边界条件，您将需要它。

至于处理狄利克雷边界条件，您可以在不首先考虑边界条件的情况下制定系统矩阵。如果给定节点的温度固定，可以这样处理：

  1. 将给定节点的右侧向量中的行设置为等于边界温度。将左侧矩阵中该行的所有列清零，并将对应于该节点的对角元素设置为等于 1。
于 2011-01-30T13:40:59.100 回答