【原】14.7 偏微分方程

14.7 偏微分方程

偏微分方程（Partial Differential Equations，简称PDEs）描述了多个变量的函数与其偏导数之间的关系。在实际应用中，很少有PDEs可以通过解析方法找到解析解，因此需要使用数值方法求解。

一类称为抛物型的偏微分方程通常以稀疏矩阵的形式得到解.

一些常用的数值方法：

1. 有限差分法（Finite Difference Method）：这是最常用的数值解PDEs的方法之一。它将PDE离散化为一个网格，并使用近似方法计算偏微分算子和导数。最常见的有限差分方法有显式差分法、隐式差分法和Crank-Nicolson方法。

2. 有限元法（Finite Element Method）：有限元法也是求解PDEs的常用方法之一。它将求解区域划分为一系列小区域，使用一组基函数（通常为多项式函数）来近似解。通过求解每个小区域内的局部问题，最终得到整个求解区域上的近似解。

3. 边界元法（Boundary Element Method）：边界元法将PDEs转化为边界积分形式。它在计算域的边界上进行离散化，并使用边界上的值和法向导数来求解方程。

4. 有限体积法（Finite Volume Method）：有限体积法将求解区域划分为一系列小体积区域，对方程进行积分。该方法基于质量守恒的原理，将PDEs转化为积分形式来求解。

这些数值方法可以应用于各种PDEs，例如抛物型、椭圆型和双曲型方程。每种方法都有其优点和适用范围，选择适当的方法要根据具体问题的性质和求解要求。

在数值求解PDEs时，还要考虑数值稳定性、收敛性和计算效率等因素。需要注意的是，对于复杂的PDEs系统，可能需要组合不同的数值方法或使用高级技术，如多重网格方法、逐步迭代求解等，以提高数值解的精度和效率。

---

在抛物线型PDE求解中，通常涉及到时间和空间两个变量。其中，时间变量是连续的，而空间变量是离散的。这种求解方法使用了时间离散化技术，如有限差分法或有限元法，将时间区间划分为若干个小时间步长。

具体举例来说，假设我们有一个一维热传导方程，它描述了材料中温度的变化随时间和空间的传播。假设材料是无限长的，其温度分布可以用函数 来表示。该方程可以写成如下的形式：

其中，α 是热传导系数。

为了求解这个方程，我们可以在空间方向上离散化 轴，将其分为一系列离散点。然后，使用有限差分法将方程转化为一个线性方程组，并使用稀疏矩阵来表示该方程组。通过求解这个线性方程组，我们可以得到温度在不同时间和空间点的数值解。

14.7.1 热传导

热沿细均匀棒的传导可以用偏微分方程来表示：

其中，表示时间 时离杆一端距离为 处的温度分布。假设沿杆长度方向没有热量损失。

解决PDE的过程中，理解符号表示法是解决问题的关键。我们建立一个矩形网格，在 和 方向上分别采用步长 和 。网格上的一个一般点具有坐标。在处用简洁的符号表示 , 记作。

然后，我们可以使用截断的泰勒级数，通过有限差分方法来近似PDE。通常情况下，方程的左侧可以使用前向差分来近似：

右侧则可以用中心差分来近似：

如果我们将方程中差分格式的右侧替换为第 行和第 行的有限差分近似的平均值，经过一定的代数运算，我们可以得到方程的以下差分格式：

其中 。这被称为Crank-Nicolson隐式方法，因为它涉及到联立方程组的解，我们将在下面看到。

为了以数值方法说明该方法，假设杆的长度为1个单位，并且它的两端与冰块接触，即边界条件为

又设初始温度(初始条件)为:

(这种情况可以通过长时间加热杆的中心，使两端与冰接触，在时间 时去除热源来实现。)这个问题在直线 处具有对称性;我们现在利用这一点来寻找解决方案。

如果取，则 ，有

将 代入该方程，可以得到以下一组方程，用于求解未知量（即经过一个时间步长 后的结果），直到杆的中点，即 ，即 。为了清晰起见，将下标 去掉：

对称性允许我们用 代替上一个方程中的 。这些方程可以写成矩阵形式

方程组左侧的矩阵 被称为三对角矩阵。解出 后，我们可以将 代入方程，然后继续求解 ，以此类推。当然，方程组可以直接在MATLAB中使用左除运算符进行求解。在下面的脚本中，方程组的一般形式被写作 。在构造矩阵 时需要小心，通常使用以下符号:

A是一个稀疏矩阵的例子.

下面的脚本实现了方程

的一般Crank-Nicolson格式，用于在 的10个时间步长内解决这个特定问题。步长由 指定，因为存在对称性。因此， 没有限制为值1，尽管此处采用了该值。脚本利用了 的稀疏性，通过使用sparse函数。

format compact
n = 5;
k = 0.01;
h = 1 / (2 * n); % 假设存在对称性
r = k / h ^ 2;
% 设置（稀疏）矩阵A
b = sparse(1:n, 1:n, 2+2*r, n, n); % b(1) .. b(n)
c = sparse(1:n-1, 2:n, -r, n, n); % c(1) .. c(n-1)
a = sparse(2:n, 1:n-1, -r, n, n); % a(2) ..
A = a + b + c;
A(n, n-1) = -2 * r; % 对称性：a(n)
full(A) %
disp(' ')
u0 = 0; % 边界条件（方程19.20）
u=2*h*[1:n]; % 初始条件（方程19.21）
u(n+1) = u(n-1); % 对称性
disp([0 u(1:n)])
for t = k*[1:10]
    g = r * ([u0 u(1:n-1)] + u(2:n+1)) + (2 - 2 * r) * u(1:n); % 方程19.19
    v = A\g'; % 方程19.24
    disp([t v'])
    u(1:n) = v;
    u(n+1) = u(n-1); % 对称性
end

注意：

• 为了保持方程的正式下标与MATLAB下标的一致性，边界值 由标量u0表示。
• 在下面的输出中，第一列是时间，后面的列是沿杆每个 时间步长的解。

0 0.2000 0.4000 0.6000 0.8000 1.0000
0.0100 0.1989 0.3956 0.5834 0.7381 0.7691
0.0200 0.1936 0.3789 0.5397 0.6461 0.6921
...
0.1000 0.0948 0.1803 0.2482 0.2918 0.3069

MATLAB中还有一些内置的PDE求解器。请参考MATLAB帮助文档：数学：微分方程：偏微分方程。