第8章数值方法8.1 数值微分8.2 数值积分8.3 常微分方程的数值解法8.4 数值求解微分方程实例

第8章数值方法

Tip

在实际工程和科学计算中，许多微积分问题无法获得解析解，或者解析表达式过于复杂。数值方法通过离散化近似，为这些问题提供了实用的解决方案。本章介绍三种基本数值方法：数值微分、数值积分和常微分方程的数值解法。

8.1 数值微分

Tip

在实际应用中，我们经常遇到需要计算函数导数的情况，但有时函数表达式复杂，或者只有离散的数据点，这时就需要使用数值微分方法。

Important

差商的基本概念

数值微分的核心思想是用差商来近似导数。差商是函数在两个点上的函数值之差与自变量之差的比值。

三种基本差商：

前向差商：
$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)}{2 h}$

$h > 0$ 是一个小的步长。

Warning

差商的几何意义

从几何角度看，差商实际上是割线斜率对切线斜率的近似：

前向差商 $(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)$ ，精度更高

Note

数值微分的计算实例

例1 $f(x) = x^2$ $x=1$ $h=0.1$

精确值 $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

可以看到中心差商给出了最精确的结果。

例2 $h$ 对精度的影响

$f(x) = \sin x$ $x=0$ $\cos 0 = 1$

$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

观察发现：

步长越小，精度越高
中心差商比前向差商精度更高
但步长过小可能引入数值舍入误差

Important

二阶导数的数值计算

二阶导数可以用以下公式近似：

f^{″} (x) \approx \frac{f (x + h) - 2 f (x) + f (x - h)}{h^{2}}

推导 $f(x+h)$ $f(x-h)$ $x$ 处泰勒展开：

f (x + h) = f (x) + h f^{'} (x) + \frac{h^{2}}{2} f^{″} (x) + \frac{h^{3}}{6} f^{‴} (x) + \dots

f (x - h) = f (x) - h f^{'} (x) + \frac{h^{2}}{2} f^{″} (x) - \frac{h^{3}}{6} f^{‴} (x) + \dots

两式相加：

f (x + h) + f (x - h) = 2 f (x) + h^{2} f^{″} (x) + O (h^{4})

整理得：

f^{″} (x) \approx \frac{f (x + h) - 2 f (x) + f (x - h)}{h^{2}}

Note

例3 $f(x) = x^3$ $x=1$ $h=0.1$

精确值 $f''(1) = 6$

数值计算：

\frac{f (1.1) - 2 f (1) + f (0.9)}{0.01} = \frac{1.331 - 2 + 0.729}{0.01} = \frac{0.06}{0.01} = 6

结果与精确值一致。

8.2 数值积分

Tip

当函数的原函数难以求出，或者只有离散的数据点时，数值积分提供了计算定积分的有效方法。数值积分的基本思想是用简单的函数（如多项式）近似被积函数，然后计算近似函数的积分。

Important

矩形法

矩形法是数值积分中最简单的方法，它将积分区间分割成若干小区间，用每个小区间左端点或右端点的函数值作为矩形的高。

左矩形法：

\int_{a}^{b} f (x) d x \approx h \sum_{i = 0}^{n - 1} f (x_{i})

$h = \frac{b-a}{n}$ $x_i = a + ih$

右矩形法：

\int_{a}^{b} f (x) d x \approx h \sum_{i = 1}^{n} f (x_{i})

几何意义：用一系列矩形面积之和近似曲边梯形的面积。

Note

例1 $\int_0^1 x^2 dx$ $n=4$

精确值 $\frac{1}{3} \approx 0.3333$

左矩形法 $h = 0.25$ $x_i = 0, 0.25, 0.5, 0.75$

$\small\begin{aligned}面积 &= 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.25 \times 0.875 = 0.21875\end{aligned}$

右矩形法 $x_i = 0.25, 0.5, 0.75, 1.0$

\begin{aligned} 面 积 & = 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.25 \times 1.875 = 0.46875 \end{aligned}

Important

梯形法

梯形法用梯形面积代替矩形面积，通常能获得更好的精度。

梯形公式：

\int_{a}^{b} f (x) d x \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) d x \approx \frac{h}{2} [f (x_{i}) + f (x_{i + 1})]

将所有小区间相加即得梯形公式。

Note

例2 $\int_0^1 x^2 dx$ $n=4$

\begin{aligned} 面 积 & = \frac{0.25}{2} [f (0) + 2 f (0.25) + 2 f (0.5) + 2 f (0.75) + f (1.0)] \\ = 0.125 \times [0 + 2 \times 0.0625 + 2 \times 0.25 + 2 \times 0.5625 + 1.0] \\ = 0.125 \times [0 + 0.125 + 0.5 + 1.125 + 1.0] = 0.125 \times 2.75 \\ = 0.34375 \end{aligned}

$0.3333$ $0.0104$ ，比矩形法有明显改进。

Important

辛普森法

辛普森法用抛物线代替直线来近似被积函数，精度更高。

辛普森公式：

\int_{a}^{b} f (x) d x \approx \frac{h}{3} [f (a) + 4 \sum_{i = 1}^{n / 2} f (x_{2 i - 1}) + 2 \sum_{i = 1}^{n / 2 - 1} f (x_{2 i}) + f (b)]

要求 $n$ 为偶数

几何意义：每两个相邻小区间上用一条抛物线近似被积函数。

Note

例3 $\int_0^1 x^2 dx$ $n=4$

\begin{aligned} 面 积 & = \frac{0.25}{3} [f (0) + 4 f (0.25) + 2 f (0.5) + 4 f (0.75) + f (1.0)] \\ = \frac{0.25}{3} \times [0 + 4 \times 0.0625 + 2 \times 0.25 + 4 \times 0.5625 + 1.0] \\ = \frac{0.25}{3} \times [0 + 0.25 + 0.5 + 2.25 + 1.0] \\ = \frac{0.25}{3} \times 4.0 = 0.3333 \end{aligned}

$x^2$ 是二次函数，辛普森法对不超过三次的多项式精确成立）。

Warning

三种数值积分方法的精度比较

$\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)$

可见辛普森法的精度最高，梯形法次之，矩形法精度最低。

Warning

数值积分的误差分析

各种数值积分方法的误差估计：

矩形法 $O(h)$
梯形法 $O(h^2)$
辛普森法 $O(h^4)$

8.3 常微分方程的数值解法

Tip

许多实际问题的数学模型是常微分方程，但大多数微分方程无法求得解析解。数值解法通过离散化方法，将连续的微分方程转化为离散的代数方程，从而获得近似解。

Important

欧拉方法

欧拉方法是最简单、最直观的数值解法，它用切线来近似积分曲线。

基本思想：从初始点出发，沿着该点的切线方向前进一小步，得到下一个近似点。

迭代公式：

y_{n + 1} = y_{n} + h f (x_{n}, y_{n})

其中：

$h$ 是步长
$f(x_n, y_n)$ $(x_n, y_n)$ 处的值
$y_n$ $x_n$ 处的近似解

几何解释：在每一步，用当前点的切线斜率代替曲线在该点的导数，沿着切线方向前进。

Note

例1 用欧拉方法求解初值问题：

\frac{d y}{d x} = y, y (0) = 1

$[0, 1]$ $h = 0.2$

精确解 $y = e^x$

数值计算：

$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

计算过程：

$y_1 = 1.0 + 0.2 \times 1.0 = 1.2$
$y_2 = 1.2 + 0.2 \times 1.2 = 1.44$
$y_3 = 1.44 + 0.2 \times 1.44 = 1.728$
依此类推...

Note

步长对精度的影响

$h$ 密切相关。一般来说，步长越小，精度越高，但计算量越大。

例2 不同步长下的比较（同一个问题）

$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

可见随着步长减小，精度逐渐提高。

Warning

改进的欧拉方法

改进的欧拉方法（又称Heun方法）通过使用预测-校正策略来提高精度。

计算步骤：

预测步：用欧拉方法计算一个初步预测值
${\bar{y}}_{n + 1} = y_{n} + h f (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})]$

几何解释：先用欧拉方法预测，然后用两端斜率的平均值来校正。

Note

例3 用改进欧拉方法求解例1的问题

第一步计算：

$\bar{y}_1 = 1.0 + 0.2 \times 1.0 = 1.2$
$y_1 = 1.0 + \frac{0.2}{2}[1.0 + 1.2] = 1.22$

$1.2214$ $1.2200$ $1.2000$ 更接近。

Note

数值解法的应用实例

应用1：人口增长模型

问题：某地区人口初始为100万，年自然增长率为3%，同时每年净迁入1万人。预测10年后的人口。

模型：

\frac{d P}{d t} = 0.03 P + 1, P (0) = 100

（单位：万人，时间单位：年）

解析解 $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电路充电过程

问题 $E = 10V$ $R = 1000\Omega$ $C = 1000\mu F$ 。求电容电压随时间的变化。

模型：

\frac{d V}{d t} = \frac{E - V}{R C}, V (0) = 0

$RC = 1000 \times 0.001 = 1$ 秒

解析解 $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$
依此类推...

Warning

数值解法的优缺点

优点：

通用性强：适用于各种类型的微分方程
实现简单：算法直观，编程容易
灵活性高：可以处理复杂的右端函数

缺点：

精度有限：欧拉方法只有一阶精度
稳定性问题：步长选择不当可能导致数值不稳定
累积误差：误差会随着计算步数增加而积累

实用建议：

先尝试欧拉方法，了解问题的基本特性
对于精度要求高的问题，使用改进欧拉方法或其他高阶方法
通过减小步长来检验结果的收敛性
对于刚性问题，需要特殊的数值方法

8.4 数值求解微分方程实例

Tip

大部分微分方程都无法解析求解, 运用计算机, 我们可以对微分方程进行数值求解.

Note

Lotka-Volterra方程

[课堂展示] 我们可以很容易的用AI帮助我们写程序, 然后画图.

Note

行星运动方程

[课堂展示] 我们可以很容易的用AI帮助我们写程序, 然后画图.

第8章 数值方法

8.1 数值微分

8.2 数值积分

8.3 常微分方程的数值解法

8.4 数值求解微分方程实例

第8章数值方法