第八章: 数值方法
在实际工程和科学计算中, 许多微积分问题无法获得解析解, 或者解析表达式过于复杂. 数值方法通过离散化近似, 为这些问题提供了实用的解决方案. 本章介绍三种基本数值方法: 数值微分, 数值积分和常微分方程的数值解法.
8.1 数值微分
在实际应用中, 我们经常遇到需要计算函数导数的情况, 但有时函数表达式复杂, 或者只有离散的数据点, 这时就需要使用数值微分方法.
数值微分的核心思想是用差商来近似导数. 差商是函数在两个点上的函数值之差与自变量之差的比值.
三种基本差商:
前向差商: f'(x) \approx \frac{f(x+h)-f(x)}{h}
后向差商: f'(x) \approx \frac{f(x)-f(x-h)}{h}
中心差商: f'(x) \approx \frac{f(x+h)-f(x-h)}{2h}
其中 h > 0 是一个小的步长.
从几何角度看, 差商实际上是割线斜率对切线斜率的近似:
前向差商: 使用点 (x, f(x)) 和 (x+h, f(x+h)) 的割线斜率
后向差商: 使用点 (x-h, f(x-h)) 和 (x, f(x)) 的割线斜率
中心差商: 使用点 (x-h, f(x-h)) 和 (x+h, f(x+h)) 的割线斜率
精度比较:
- 前向和后向差商的误差约为 O(h)
- 中心差商的误差约为 O(h^2), 精度更高
计算 f(x) = x^2 在 x=1 处的导数近似值，取步长 h=0.1，分别用前向差商、后向差商和中心差商，并与精确值 f'(1) = 2 比较。
三种差商公式分别为：前向 (f(x+h)-f(x))/h，后向 (f(x)-f(x-h))/h，中心 (f(x+h)-f(x-h))/(2h)。
精确值：f'(1) = 2。
前向差商：
\frac{f(1.1)-f(1)}{0.1} = \frac{1.21-1}{0.1} = 2.1.
后向差商：
\frac{f(1)-f(0.9)}{0.1} = \frac{1-0.81}{0.1} = 1.9.
中心差商：
\frac{f(1.1)-f(0.9)}{0.2} = \frac{1.21-0.81}{0.2} = 2.0.
中心差商给出了最精确的结果，与精确值完全一致，因为 f(x)=x^2 是二次函数，中心差商对二次函数精确成立。
对于 f(x) = \sin x，在 x=0 处计算一阶导数（精确值 \cos 0 = 1），分析步长 h 对前向差商和中心差商精度的影响。
分别取 h = 0.1, 0.01, 0.001，计算两种差商并与精确值对比，观察误差随步长减小的变化趋势。
| 步长 h | 前向差商 | 绝对误差 | 中心差商 | 绝对误差 |
|---|---|---|---|---|
| 0.1 | 0.99833 | 0.00167 | 0.99998 | 0.00002 |
| 0.01 | 0.99998 | 0.00002 | 1.00000 | 0.00000 |
| 0.001 | 1.00000 | 0.00000 | 1.00000 | 0.00000 |
结论： - 步长越小，精度越高。 - 中心差商比前向差商精度更高（误差约为 O(h^2) vs O(h)）。 - 但步长过小可能引入数值舍入误差。
二阶导数可以用以下公式近似: f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}
推导: 将 f(x+h) 和 f(x-h) 在 x 处泰勒展开: f(x+h) = f(x) + hf'(x) + \frac{h^2}{2}f''(x) + \frac{h^3}{6}f'''(x) + \cdots f(x-h) = f(x) - hf'(x) + \frac{h^2}{2}f''(x) - \frac{h^3}{6}f'''(x) + \cdots
两式相加: f(x+h) + f(x-h) = 2f(x) + h^2 f''(x) + O(h^4)
整理得: f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}
计算 f(x) = x^3 在 x=1 处的二阶导数，取 h=0.1，并与精确值 f''(1) = 6 比较。
使用二阶导数数值公式 f''(x) \approx [f(x+h) - 2f(x) + f(x-h)]/h^2。
\frac{f(1.1) - 2f(1) + f(0.9)}{h^2} = \frac{1.331 - 2 + 0.729}{0.01} = \frac{0.06}{0.01} = 6.
结果与精确值 f''(1) = 6 完全一致。（f(x)=x^3 是三次函数，该公式对三次函数精确成立。）
8.2 数值积分
当函数的原函数难以求出, 或者只有离散的数据点时, 数值积分提供了计算定积分的有效方法. 数值积分的基本思想是用简单的函数 (如多项式) 近似被积函数, 然后计算近似函数的积分.
矩形法是数值积分中最简单的方法, 它将积分区间分割成若干小区间, 用每个小区间左端点或右端点的函数值作为矩形的高.
左矩形法: \int_a^b f(x)dx \approx h\sum_{i=0}^{n-1} f(x_i) 其中 h = \frac{b-a}{n}, x_i = a + ih
右矩形法: \int_a^b f(x)dx \approx h\sum_{i=1}^n f(x_i)
几何意义: 用一系列矩形面积之和近似曲边梯形的面积.
用左矩形法和右矩形法计算 \displaystyle\int_0^1 x^2\,dx，取 n=4，并与精确值 \dfrac{1}{3} \approx 0.3333 比较。
h = 0.25，左矩形法取各区间左端点，右矩形法取各区间右端点，乘以步长后求和。
左矩形法（x_i = 0, 0.25, 0.5, 0.75）：
0.25 \times [f(0) + f(0.25) + f(0.5) + f(0.75)] = 0.25 \times [0 + 0.0625 + 0.25 + 0.5625] = 0.21875.
右矩形法（x_i = 0.25, 0.5, 0.75, 1.0）：
0.25 \times [f(0.25) + f(0.5) + f(0.75) + f(1.0)] = 0.25 \times [0.0625 + 0.25 + 0.5625 + 1.0] = 0.46875.
左矩形法偏小，右矩形法偏大，两者均与精确值有较大误差。
梯形法用梯形面积代替矩形面积, 通常能获得更好的精度.
梯形公式: \int_a^b f(x)dx \approx \frac{h}{2}[f(a) + 2\sum_{i=1}^{n-1} f(x_i) + f(b)]
几何意义: 用一系列梯形面积之和近似曲边梯形的面积.
推导: 每个小区间 [x_i, x_{i+1}] 上的积分用梯形面积近似: \int_{x_i}^{x_{i+1}} f(x)dx \approx \frac{h}{2}[f(x_i) + f(x_{i+1})]
将所有小区间相加即得梯形公式.
用梯形法计算 \displaystyle\int_0^1 x^2\,dx，取 n=4，并与精确值 \dfrac{1}{3} 比较。
梯形公式：\dfrac{h}{2}[f(a) + 2\sum_{i=1}^{n-1}f(x_i) + f(b)]，h=0.25。
\frac{0.25}{2}[f(0) + 2f(0.25) + 2f(0.5) + 2f(0.75) + f(1.0)]
= 0.125 \times [0 + 2(0.0625) + 2(0.25) + 2(0.5625) + 1.0] = 0.125 \times 2.75 = 0.34375.
与精确值 0.3333 相比，误差约 0.0104，比矩形法有明显改进。
辛普森法用抛物线代替直线来近似被积函数, 精度更高.
辛普森公式: \int_a^b f(x)dx \approx \frac{h}{3}[f(a) + 4\sum_{i=1}^{n/2} f(x_{2i-1}) + 2\sum_{i=1}^{n/2-1} f(x_{2i}) + f(b)]
要求: 区间等分且 n 为偶数
几何意义: 每两个相邻小区间上用一条抛物线近似被积函数.
用辛普森法计算 \displaystyle\int_0^1 x^2\,dx，取 n=4，并与精确值 \dfrac{1}{3} 比较。
辛普森公式：\dfrac{h}{3}[f(a) + 4f(x_1) + 2f(x_2) + 4f(x_3) + f(b)]（n=4，h=0.25）。
\frac{0.25}{3}[f(0) + 4f(0.25) + 2f(0.5) + 4f(0.75) + f(1.0)]
= \frac{0.25}{3}[0 + 4(0.0625) + 2(0.25) + 4(0.5625) + 1.0] = \frac{0.25}{3} \times 4.0 = 0.3333.
与精确值完全一致，因为 x^2 是二次多项式，辛普森法对不超过三次的多项式精确成立。
对于 \int_0^1 e^x dx = e - 1 \approx 1.71828, 不同方法的精度比较:
| 方法 | n=4 时的近似值 | 绝对误差 | 收敛阶 |
|---|---|---|---|
| 左矩形法 | 1.49068 | 0.22760 | O(h) |
| 右矩形法 | 1.93637 | 0.21809 | O(h) |
| 梯形法 | 1.71339 | 0.00489 | O(h^2) |
| 辛普森法 | 1.71832 | 0.00004 | O(h^4) |
可见辛普森法的精度最高, 梯形法次之, 矩形法精度最低.
各种数值积分方法的误差估计:
- 矩形法: 误差约 O(h)
- 梯形法: 误差约 O(h^2)
- 辛普森法: 误差约 O(h^4)
8.3 常微分方程的数值解法
许多实际问题的数学模型是常微分方程, 但大多数微分方程无法求得解析解. 数值解法通过离散化方法, 将连续的微分方程转化为离散的代数方程, 从而获得近似解.
欧拉方法是最简单, 最直观的数值解法, 它用切线来近似积分曲线.
基本思想: 从初始点出发, 沿着该点的切线方向前进一小步, 得到下一个近似点.
迭代公式: y_{n+1} = y_n + hf(x_n, y_n)
其中: - h 是步长 - f(x_n, y_n) 是微分方程右端函数在 (x_n, y_n) 处的值 - y_n 是 x_n 处的近似解
几何解释: 在每一步, 用当前点的切线斜率代替曲线在该点的导数, 沿着切线方向前进.
用欧拉方法求解初值问题
\frac{dy}{dx} = y, \quad y(0) = 1
在区间 [0, 1] 上，取步长 h = 0.2，精确解为 y = e^x。
欧拉迭代公式：y_{n+1} = y_n + h \cdot f(x_n, y_n)，其中 f(x, y) = y。
迭代过程：
y_1 = 1.0 + 0.2 \times 1.0 = 1.2,\quad y_2 = 1.2 + 0.2 \times 1.2 = 1.44,
y_3 = 1.44 + 0.2 \times 1.44 = 1.728,\quad y_4 = 1.728 + 0.2 \times 1.728 = 2.0736,
y_5 = 2.0736 + 0.2 \times 2.0736 = 2.4883.
| x_n | 欧拉近似 y_n | 精确值 e^{x_n} | 绝对误差 |
|---|---|---|---|
| 0.0 | 1.0000 | 1.0000 | 0.0000 |
| 0.2 | 1.2000 | 1.2214 | 0.0214 |
| 0.4 | 1.4400 | 1.4918 | 0.0518 |
| 0.6 | 1.7280 | 1.8221 | 0.0941 |
| 0.8 | 2.0736 | 2.2255 | 0.1519 |
| 1.0 | 2.4883 | 2.7183 | 0.2300 |
误差随 x 增大而积累，步长越小精度越高。
对初值问题 \dfrac{dy}{dx} = y，y(0) = 1，用欧拉方法分别取步长 h = 0.5, 0.2, 0.1, 0.05, 0.01，比较 y(1) 的近似值与精确值 e \approx 2.7183 的误差。
对每个步长反复应用欧拉迭代公式至 x=1，记录最终近似值并计算误差。
| 步长 h | y(1) 近似值 | 绝对误差 |
|---|---|---|
| 0.5 | 2.2500 | 0.4683 |
| 0.2 | 2.4883 | 0.2300 |
| 0.1 | 2.5937 | 0.1246 |
| 0.05 | 2.6533 | 0.0650 |
| 0.01 | 2.7048 | 0.0135 |
随着步长减小，精度逐渐提高，误差约为 O(h)。
改进的欧拉方法 (又称Heun方法) 通过使用预测-校正策略来提高精度.
计算步骤: 1. 预测步: 用欧拉方法计算一个初步预测值 \bar{y}_{n+1} = y_n + hf(x_n, y_n)
- 校正步: 用预测点的斜率来改进 y_{n+1} = y_n + \frac{h}{2}[f(x_n, y_n) + f(x_{n+1}, \bar{y}_{n+1})]
几何解释: 先用欧拉方法预测, 然后用两端斜率的平均值来校正.
用改进欧拉方法（Heun方法）求解初值问题 \dfrac{dy}{dx} = y，y(0) = 1，取步长 h = 0.2，计算第一步的近似值，并与普通欧拉方法和精确值 e^{0.2} \approx 1.2214 比较。
改进欧拉法：先预测 \bar{y}_1 = y_0 + h f(x_0, y_0)，再校正 y_1 = y_0 + \dfrac{h}{2}[f(x_0,y_0) + f(x_1,\bar{y}_1)]。
预测步：
\bar{y}_1 = 1.0 + 0.2 \times 1.0 = 1.2.
校正步：
y_1 = 1.0 + \frac{0.2}{2}[f(0, 1.0) + f(0.2, 1.2)] = 1.0 + 0.1[1.0 + 1.2] = 1.22.
改进欧拉法结果 1.22 比普通欧拉法的 1.20 更接近精确值 1.2214，精度有所提升。
应用1：人口增长模型
某地区人口初始为100万，年自然增长率为3%，同时每年净迁入1万人。预测10年后的人口。
模型： \frac{dP}{dt} = 0.03P + 1, \quad P(0) = 100 （单位：万人，时间单位：年）
用欧拉方法（h = 1 年）计算前5年的人口，并与解析解 P(t) = \frac{100}{3}(e^{0.03t} - 1) + 100e^{0.03t} 比较。
应用2：RC电路充电过程
RC电路中，电容初始电压为0，电源电压 E = 10V，电阻 R = 1000\Omega，电容 C = 1000\mu F。求电容电压随时间的变化。
模型： \frac{dV}{dt} = \frac{E - V}{RC}, \quad V(0) = 0
其中 RC = 1000 \times 0.001 = 1 秒。
用欧拉方法（h = 0.2 秒）计算前5步，并与解析解 V(t) = 10(1 - e^{-t}) 比较。
欧拉方法迭代公式：y_{n+1} = y_n + h \cdot f(t_n, y_n)。
应用1：人口增长模型
解析解：P(t) = \frac{100}{3}(e^{0.03t} - 1) + 100e^{0.03t}
数值解（欧拉方法，h = 1 年）：
| 年份 t | 人口 P(t)（数值） | 人口 P(t)（精确） |
|---|---|---|
| 0 | 100.00 | 100.00 |
| 1 | 104.00 | 104.08 |
| 2 | 108.12 | 108.37 |
| 3 | 112.36 | 112.78 |
| 4 | 116.73 | 117.33 |
| 5 | 121.23 | 122.02 |
| 10 | 147.85 | 149.18 |
计算过程： - \small P_1 = 100 + 1 \times (0.03 \times 100 + 1) = 104 - \small P_2 = 104 + 1 \times (0.03 \times 104 + 1) = 108.12 - 依此类推…
应用2：RC电路充电过程
解析解：V(t) = 10(1 - e^{-t})
数值解（欧拉方法，h = 0.2 秒）：
| 时间 t (s) | 电压 V(t)（数值） | 电压 V(t)（精确） |
|---|---|---|
| 0.0 | 0.00 | 0.00 |
| 0.2 | 2.00 | 1.81 |
| 0.4 | 3.60 | 3.30 |
| 0.6 | 4.88 | 4.51 |
| 0.8 | 5.90 | 5.51 |
| 1.0 | 6.72 | 6.32 |
计算过程： - V_1 = 0 + 0.2 \times (10 - 0)/1 = 2.0 - V_2 = 2.0 + 0.2 \times (10 - 2.0)/1 = 3.6 - 依此类推…
优点: 1. 通用性强: 适用于各种类型的微分方程 2. 实现简单: 算法直观, 编程容易 3. 灵活性高: 可以处理复杂的右端函数
缺点: 1. 精度有限: 欧拉方法只有一阶精度 2. 稳定性问题: 步长选择不当可能导致数值不稳定 3. 累积误差: 误差会随着计算步数增加而积累
实用建议: 1. 先尝试欧拉方法, 了解问题的基本特性 2. 对于精度要求高的问题, 使用改进欧拉方法或其他高阶方法 3. 通过减小步长来检验结果的收敛性 4. 对于刚性问题, 需要特殊的数值方法
8.4 数值求解微分方程实例
大部分微分方程都无法解析求解, 运用计算机, 我们可以对微分方程进行数值求解.
[课堂展示] 我们可以很容易的用AI帮助我们写程序, 然后画图.
[课堂展示] 我们可以很容易的用AI帮助我们写程序, 然后画图.