下面的代码求解一维热方程,该方程代表一根杆,其末端保持在零温度,初始条件为 10*np.sin(np.pi*x)。
Dirichlet 边界条件(两端为零温度)如何应用于此计算?有人告诉我矩阵 A 的上下两行包含两个非零元素,而缺少的第三个元素 是 Dirichlet 条件。但我不明白这种情况是通过哪种机制影响计算的。 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)
最佳答案
让我们看一个简单的例子。我们假设 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
但是,由于齐次 Dirichlet 边界条件,所以我们可以简单地从等式中省略它,并且 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 /
矩阵中缺失的条目确实对应于边界条件为零的事实。
关于python - 边界条件在热方程和Crank-Nicholson有限差分解中的应用,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/4843034/