微分方程一维抛物热传导方程数值解法全解析-智慧文博士

微分方程一维抛物热传导方程向前向后欧拉C-N格式二阶BDF格式MATLAB源码显式欧拉，隐式欧拉，梯形公式，改进欧拉五点差分，九点差分差分格式，紧差分格式直拍，只有pdf版方法说明 word版公式纯手打数值例子有数据图解分析含源码和流程图

在科学与工程领域，一维抛物热传导方程是描述热传递等诸多现象的重要模型。今天咱就来唠唠求解它的几种经典数值方法，从显式欧拉、隐式欧拉这些基础的，到梯形公式、改进欧拉，还有差分格式里的五点差分、九点差分以及紧差分格式，最后咱还会给出MATLAB源码，以及带数据图解分析的数值例子。

1. 显式欧拉法

显式欧拉法是求解微分方程的基础方法。对于一维抛物热传导方程：

\[ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \]

离散化后，显式欧拉的时间推进公式为：

\[ uj^{n + 1} = uj^n + \frac{\alpha \Delta t}{\Delta x^2} (u{j - 1}^n - 2uj^n + u_{j + 1}^n) \]

这里 \( u_j^n \) 表示在 \( n \) 时刻， \( j \) 空间位置的温度值。

MATLAB 代码示例：

% 显式欧拉法求解一维热传导方程 L = 1; % 区域长度 T = 1; % 总时间 Nx = 100; % 空间节点数 Nt = 1000; % 时间节点数 dx = L / (Nx - 1); dt = T / Nt; alpha = 1; % 热扩散系数 u = zeros(Nx, Nt + 1); x = linspace(0, L, Nx); t = linspace(0, T, Nt + 1); % 初始条件 u(:, 1) = sin(pi * x); for n = 1:Nt for j = 2:Nx - 1 u(j, n + 1) = u(j, n) + alpha * dt / dx^2 * (u(j - 1, n) - 2 * u(j, n) + u(j + 1, n)); end % 边界条件 u(1, n + 1) = 0; u(Nx, n + 1) = 0; end % 绘图 figure; for n = 1:10:Nt + 1 plot(x, u(:, n)); hold on; end xlabel('x'); ylabel('u(x, t)'); title('显式欧拉法求解一维热传导方程'); hold off;

代码分析：首先定义了区域长度、总时间、空间和时间节点数等参数。通过循环来实现时间推进，在每个时间步内，根据显式欧拉公式更新内部节点的值，同时考虑边界条件。最后绘图展示不同时刻温度分布。

2. 隐式欧拉法

隐式欧拉法在稳定性上比显式欧拉更具优势。其离散化公式为：

\[ uj^{n + 1} - \frac{\alpha \Delta t}{\Delta x^2} (u{j - 1}^{n + 1} - 2uj^{n + 1} + u{j + 1}^{n + 1}) = u_j^n \]

这是一个关于 \( u_j^{n + 1} \) 的线性方程组，需要求解。

MATLAB 代码如下：

% 隐式欧拉法求解一维热传导方程 L = 1; T = 1; Nx = 100; Nt = 1000; dx = L / (Nx - 1); dt = T / Nt; alpha = 1; u = zeros(Nx, Nt + 1); x = linspace(0, L, Nx); t = linspace(0, T, Nt + 1); u(:, 1) = sin(pi * x); % 构建系数矩阵 A = spdiags([ones(Nx - 1, 1), -2 * ones(Nx, 1), ones(Nx - 1, 1)], [-1, 0, 1], Nx, Nx); A(1, 1) = 1; A(Nx, Nx) = 1; A = A - (alpha * dt / dx^2) * A; for n = 1:Nt b = u(:, n); b(1) = 0; b(Nx) = 0; u(:, n + 1) = A \ b; end figure; for n = 1:10:Nt + 1 plot(x, u(:, n)); hold on; end xlabel('x'); ylabel('u(x, t)'); title('隐式欧拉法求解一维热传导方程'); hold off;

代码分析：这里先构建了系数矩阵 \( A \)，它与隐式欧拉的离散化方程相关。在每个时间步，通过求解线性方程组 \( A u^{n + 1} = b \) 来得到下一个时间步的温度分布。

3. 梯形公式（C - N 格式）

梯形公式又称 Crank - Nicolson 格式，它是一种隐式的二阶精度方法。离散化公式为：

\[ uj^{n + 1} - \frac{\alpha \Delta t}{2\Delta x^2} (u{j - 1}^{n + 1} - 2uj^{n + 1} + u{j + 1}^{n + 1}) = uj^n + \frac{\alpha \Delta t}{2\Delta x^2} (u{j - 1}^n - 2uj^n + u{j + 1}^n) \]

MATLAB 代码：

% Crank - Nicolson格式求解一维热传导方程 L = 1; T = 1; Nx = 100; Nt = 1000; dx = L / (Nx - 1); dt = T / Nt; alpha = 1; u = zeros(Nx, Nt + 1); x = linspace(0, L, Nx); t = linspace(0, T, Nt + 1); u(:, 1) = sin(pi * x); % 构建系数矩阵 A = spdiags([ones(Nx - 1, 1), -2 * ones(Nx, 1), ones(Nx - 1, 1)], [-1, 0, 1], Nx, Nx); A(1, 1) = 1; A(Nx, Nx) = 1; A = A - (alpha * dt / (2 * dx^2)) * A; for n = 1:Nt b = u(:, n) + (alpha * dt / (2 * dx^2)) * (A * u(:, n)); b(1) = 0; b(Nx) = 0; u(:, n + 1) = A \ b; end figure; for n = 1:10:Nt + 1 plot(x, u(:, n)); hold on; end xlabel('x'); ylabel('u(x, t)'); title('Crank - Nicolson格式求解一维热传导方程'); hold off;

代码分析：同样是构建系数矩阵 \( A \)，但这里矩阵与时间步的右侧项 \( b \) 构建都与梯形公式相关。通过求解线性方程组来更新温度分布。

4. 改进欧拉法

改进欧拉法对显式欧拉进行了改进，提高了精度。它先通过显式欧拉预估，再进行校正。

五点差分与九点差分

五点差分是常用的二阶精度差分格式，对于二阶导数 \(\frac{\partial^2 u}{\partial x^2} \) 的五点差分近似为：

\[ \frac{\partial^2 u}{\partial x^2} \approx \frac{-u{j + 2} + 16u{j + 1} - 30uj + 16u{j - 1} - u_{j - 2}}{12\Delta x^2} \]

九点差分格式则在五点的基础上利用更多节点来提高精度，公式相对复杂些。

紧差分格式

紧差分格式利用相邻节点的线性组合来逼近导数，能在较少节点下达到较高精度。

数值例子与数据图解分析

咱以一个具体的例子来看，比如初始条件 \( u(x, 0) = \sin(\pi x) \)，边界条件 \( u(0, t) = u(1, t) = 0 \)。通过上述不同方法求解后，我们可以绘制不同时刻温度分布曲线。从图中可以直观看到显式欧拉法在时间步较大时可能出现不稳定，而隐式方法和 C - N 格式更加稳定。不同差分格式在精度上也有差异，紧差分格式在某些情况下能以较少节点达到较高精度。