经典数值计算

插值问题

物理问题 $\rightarrow$ 离散化 $\rightarrow$ 计算机处理

离散的物理数据，但不知道实际函数$\rightarrow$ 构造近似函数作为实际函数$f(x)$的逼近

---

插值问题的数学表述： $f(x)$是定义在区间$[a, b]$的函数，$(x_0, y_0),\ (x_1, y_1),\ (x_2, y_2)...,\ (x_n, y_n)$是该区间上$n+1$个不同的点。我们需要构造一个便于计算的函数$P(x)$，s.t.

$$P(x_i) = y_i$$

则称$P(x)$是$f(x)$的插值函数，$[a, b]$是插值区间，$(x_i, y_i)$是插值节点。同时在其他$x \neq x_i$点上，$P(x)$可以作为$f(x)$的近似。

多项式插值

设$P_n(x) = a_0 + a_1x + a_2x^2 + ... + a_nx^n,\;\; a_0, a_1, ..., a_n \in \mathbb{R}$

---

定理：对于$n+1$个互不相同的插值节点，满足插值条件的$n$次多项式插值函数存在且唯一。

---

多项式插值方法一 —— 拉格朗日插值法

引入：线性插值

给定两个点$(x_0, y_0)$，$(x_1, y_1)$求解

$$\begin{cases} a_0+ a_1x_0 = y_0\\ a_0+ a_1x_1 = y_1 \end{cases}$$

可得

$$P_1(x) = a_0+ a_1x = \frac{x-x_1}{x_0 - x_1}y_0 + \frac{x-x_0}{x_1 - x_0}y_1$$

推广：拉格朗日插值

$$P_n(x) = \sum_{k=0}^n l_k(x) \cdot y_k, \quad l_k(x) = \prod_{j \neq k} \frac{x - x_j}{x_k - x_j}$$

其中$l_k(x)$被称作基函数。

拉格朗日插值基函数满足以下性质：

• 归一性：$l_k(x_i) = \begin{cases} 1 \quad i = k\\ 0 \quad o/w \end{cases}$
• 多项式次数：每个$l_k$是$n$次多项式。

举例理解

上面的线性插值已经是一个例子。现在看$n = 2$：

对于$(x_0,x_1,x_2)$，

$$\begin{aligned} l_0(x) &= \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} \\ l_1(x) &= \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} \\ l_2(x) &= \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} \end{aligned}$$

线性插值误差

---

定理：设 $f(x)$ 一阶连续可导，且 $f''(x)$ 在区间 $(a, b)$ 上存在，则对任意 $x \in [a, b]$，存在 $\xi \in [a, b]$ 使得：

$$R(x) = f(x) - P_1(x) = \frac{f''(\xi)}{2}(x-x_0)(x-x_1)$$

---

证明过程：

定义误差函数：

$$R(x) = f(x) - P_1(x) = k(x)(x-x_0)(x-x_1)$$

构造辅助函数：

$$\phi(t) = f(t) - P_1(t) - k(x)(t-x_0)(t-x_1)$$

$\phi(t)$ 在 $x_0,x,x_1$ 处有零点，且$\phi'(t)$ 在 $(\xi_1,\xi_2)$ 有零点

$\therefore$ 应用Rolle定理：$\phi''(\xi) = 0$ $\Rightarrow$ $f''(\xi) = 2k(x)$

一般形式的误差估计： 对于 $n$ 次插值多项式 $P_n(x)$：

$$R_n(x) = f(x) - P_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)$$

其中：

$$\omega_{n+1}(x) = \prod_{i=0}^n (x-x_i)$$

$\therefore$

1. 误差与高阶导数 $f^{(n+1)}$ 成正比
2. 节点越密集 ($(x-x_i)$ 越小)，误差越小
3. 增加节点数 $n$ 可提高精度

很多情况下，$f^{(n+1)}$ 未知，可采用以下方法估计：

使用 $n+2$ 个数据点进行两次插值运算：前 $n+1$ 点构造 $P_n(x)$；后 $n+1$ 点构造 $\hat{P}_n(x)$。整理可得误差关系式：

$$f(x) - P_n(x) \approx \frac{x-x_0}{x_0-x_{n+1}}(P_n(x)-\hat{P}_n(x))$$

拉格朗日插值法的优点/缺点：

• 优点：
  • 显示公式：无需解线性方程组，直接构造多项式。
  • 理论直观
• 缺点：
  • 计算效率低：新增节点需重新计算所有基函数。
  • 数值不稳定：高次插值时可能因舍入误差放大失效。

多项式插值方法二——牛顿插值法

引入

避开了拉格朗日插值法“新增节点需重新计算所有基函数”这一问题。

由数学知识可知任何一个多项式都可以表示为：

$$N_n(x) = a_0 + \sum_{k=1}^n a_k \prod_{i=0}^{k-1} (x - x_i)$$

在牛顿形式下，我们可以认为此时$a_k$对应的基函数是$k$的多项式$\prod_{i=0}^{k-1} (x - x_i)$，它们相互独立。因此有更高阶数加入的时候只需计算相应新的多项式系数即可。

计算步骤

将插值条件$P(x_i) = y_i$写成线性方程组的形式：

$\begin{cases} y_0 = a_0 = f(x_0)\\ y_1 = y_0 +a_1 (x_1-x_0) = f(x_1)\\ y_2 = y_0 + a_1 (x_2-x_0) + a_2(x_2-x_0)(x_2-x_1) = f(x_2)\\ ...\\ y_n = y_0 + a_1 (x_n-x_0) + a_2(x_n-x_0)(x_n-x_1) + ... + a_n(x_n-x_0)(x_n-x_1)...(x_n-x_{n-1}) = f(x_n) \end{cases}$
引入差商的概念：
由$(y_0, y_1)$我们引入一阶差商

$$F(x_0, x_1) = \frac{f(x_1) - f(x_0)}{x_1 - x_0}$$

类似可得更高阶的差商公式：

$$F(x_2, x_1) = \frac{f(x_2) - f(x_1)}{x_2 - x_1}\\ F(x_0, x_1,x_2) = \frac{F(x_2, x_1)-F(x_0, x_1)}{x_2-x_0}\\ ...\\ F(x_0,...,x_n) = \frac{F(x_1,...,x_n)-F(x_0,...,x_{n-1})}{x_n-x_0}$$

根据以上公式我们可以改写差商形式的牛顿插值多项式：

$$P_n(x) = F(x_0) + \sum_{k=1}^n F(x_0,...,x_k) \cdot \prod_{i=0}^{k-1}(x-x_i)$$

示例

$$\begin{aligned} P_3(x) = & F(x_0)+ F(x_0,x_1)(x-x_0) \\ & + F(x_0,x_1,x_2)(x-x_0)(x-x_1) \\ & + F(x_0,x_1,x_2,x_3)(x-x_0)(x-x_1)(x-x_2) \end{aligned}$$

牛顿差商的性质

1. 线性性：
   若$h(x) = \alpha f(x) + \beta g(x)$，则：

$$H(x_0,...,x_n) = \alpha F(x_0,...,x_n) + \beta G(x_0,...,x_n)$$

2. 对称性
   差商值与节点排列顺序无关：

$$F(x_0,x_1,...,x_n) = F(x_{\sigma(0)},x_{\sigma(1)},...,x_{\sigma(n)})$$

其中$\sigma$是任意排列

3. 差商和导数关系：
   若$f(x)\in C^n[a,b]$，则存在$\xi\in(a,b)$使得：

4. 差商表可加性：
   新增节点$x_{n+1}$时：

$$F(x_0,...,x_{n+1}) = \frac{F(x_1,...,x_{n+1}) - F(x_0,...,x_n)}{x_{n+1}-x_0}$$

$$F(x_0,...,x_n) = \frac{f^{(n)}(\xi)}{n!}$$

荣格(Runge)现象

指在高次多项式插值中，随着插值节点数量的增加，插值结果在区间端点附近出现剧烈振荡的现象。这种现象由德国数学家Carl Runge在研究多项式插值时首次发现。

解决方法：多次线性插值

将插值区间划分为若干子区间，在每个子区间上用线性多项式进行插值的方法。给定节点$(x_0,y_0),(x_1,y_1),...,(x_n,y_n)$，我们构造的插值函数应该满足：

• $S(x_i) = y_i$ （通过所有节点）
• $S(x)$在每个$[x_i,x_{i+1}]$上是线性函数
• 整体函数连续（$C^0$连续）

具体实现步骤：

1. 找到待插值点$x$及其所在的区间$[x_k,x_{k+1}]$
2. 计算：$ P_1(x) = \frac{x-x_{k+1}}{x_k - x_{k+1}}y_k + \frac{x-x_k}{x_{k+1} - x_k}y_{k+1} $

缺陷：$C^1$不连续

样条插值

样条插值（Spline Interpolation）是一种分段低次多项式插值方法，通过强制拼接点处的光滑条件来保证整体曲线的平滑性。

插值函数需要满足：

• $S_i(x_i)=y_i$, $S_i(x_{i+1})=y_{i+1}$（插值条件）
• 连续性：
  • $C^0$: $S_{i-1}(x_i) = S_i(x_i)$
  • $C^1$: $S'_{i-1}(x_i) = S'_i(x_i)$
  • $C^2$: $S''_{i-1}(x_i) = S''_i(x_i)$

最常用：三次样条插值函数，即在每段$[x_i,x_{i+1}]$上构造：

$$S_i(x) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3$$

对$n+1$个插值点，总共有$n$个区间，那么总参数有$4n$个。

约束条件：

• 插值条件：提供 $n+1$ 个
• 连续性：提供 $3(n-1)$ 个
• 边界条件：提供n个

对于$x \in [x_k, x_{k+1})$：

$$\begin{aligned} S_k(x) &= a_k + b_k h + c_k h^2 + d_k h^3 \quad (h = x-x_k) \\ S'_k(x) &= b_k + 2c_k h + 3d_k h^2 \\ S''_k(x) &= 2c_k + 6d_k h \end{aligned}$$

应满足条件：

• $S_k(x_k) = y_k \Rightarrow a_k = y_k$
• $ S_k(x_{k+1}) = y_{k+1} \Rightarrow a_k + h_k b_k + h_k^2 c_k + h_k^3 d_k = y_{k+1} $ （其中 $h_k = x_{k+1}-x_k$）
• $ b_k + 2h_k c_k + 3h_k^2 d_k = b_{k+1} $
• $ 2c_k + 6h_k d_k = 2c_{k+1} $

求解该方程组：
令 $m_k = 2c_k$，则：

1. 由二阶导数连续条件：

   $$d_k = \frac{m_{k+1}-m_k}{6h_k}$$

2. 代入插值方程得：

   $$b_k = \frac{y_{k+1}-y_k}{h_k} - \frac{h_k}{2}m_k - \frac{h_k}{6}(m_{k+1}-m_k)$$

3. 将系数代入一阶导数连续条件，得到：

$$h_k m_k + 2(h_k+h_{k+1})m_{k+1} + h_{k+1}m_{k+2} = 6\left(\frac{y_{k+2}-y_{k+1}}{h_{k+1}} - \frac{y_{k+1}-y_k}{h_k}\right)$$

4. 将上两个边界条件的约束 $\Rightarrow$ 建立起线性方程组即可求解。

数据拟合

数据拟合是通过数学模型近似描述离散数据点的方法，与插值的区别在于：

• 插值：强制通过所有数据点
• 拟合：允许存在偏差，追求整体趋势匹配

最小二乘法拟合

数学原理：最小化残差平方和：

$$\min \sum_{i=1}^n (y_i - f(x_i))^2$$

线性拟合：

考虑$f(x) = a + bx$。
残差平方和：$$Q(a,b) = \sum_{i=1}^n (f(x_i) - y_i)^2 = \sum_{i=1}^n (a + b x_i - y_i)^2$$
$Q(a, b)$最小，需满足：

$$\begin{cases} \frac{\partial Q}{\partial a} = 2\sum_{i=1}^n (a + b x_i - y_i) = 0\\ \frac{\partial Q}{\partial b} = 2\sum_{i=1}^n (a + b x_i - y_i)x_i = 0 \end{cases}$$

将方程组整理为矩阵形式：

$$\begin{bmatrix} n & \sum x_i \\ \sum x_i & \sum x_i^2 \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \end{bmatrix}$$

求解$(a, b)$即可。

多项式拟合

k次多项式模型：

$$f(x) = \sum_{j=0}^k a_j x^j$$

残差平方和：

$$\begin{cases} a_0 \sum 1 + a_1 \sum x_i + \cdots + a_m \sum x_i^m = \sum y_i \\ a_0 \sum x_i + a_1 \sum x_i^2 + \cdots + a_m \sum x_i^{m+1} = \sum x_i y_i \\ \vdots \\ a_0 \sum x_i^m + a_1 \sum x_i^{m+1} + \cdots + a_m \sum x_i^{2m} = \sum x_i^m y_i \end{cases}$$

整理成矩阵形式：

$$\begin{bmatrix} n & \sum x_i & \cdots & \sum x_i^m \\ \sum x_i & \sum x_i^2 & \cdots & \sum x_i^{m+1} \\ \vdots & \vdots & \ddots & \vdots \\ \sum x_i^m & \sum x_i^{m+1} & \cdots & \sum x_i^{2m} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_m \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \\ \vdots \\ \sum x_i^m y_i \end{bmatrix}$$

求解即可。

非线性拟合

模型类型	变换方法	线性形式
指数模型 $y=ae^{bx}$	取对数	$\ln y = \ln a + bx$
幂函数 $y=ax^b$	取对数	$\ln y = \ln a + b\ln x$
对数模型 $y=a\ln x + b$	-	直接拟合
倒数模型 $y = \frac{x}{a+bx}$	取倒数	$\frac{1}{y} = \frac{a}{x} + b$

容易出现的问题：欠拟合&过拟合。

数值微分

用离散点进行数值微分：

假设数据点均匀分布，$x_i = x_0 + i \cdot h$，在$x_{i-1}$、$x_i$、$x_{i+1}$附近Taylor展开：

$$f_{i+1} = f_i + h f'_i + \frac{h^2}{2} f''_i + \frac{h^3}{6} f'''_i + O(h^4)$$

$$f_{i-1} = f_i - h f'_i + \frac{h^2}{2} f''_i - \frac{h^3}{6} f'''_i + O(h^4)$$

就可以得到以下公式：

类型	公式	误差阶
前向差分	$f'(x) \approx \frac{f(x+h)-f(x)}{h}$	$O(h)$
后向差分	$f'(x) \approx \frac{f(x)-f(x-h)}{h}$	$O(h)$
中心差分	$f'(x) \approx \frac{f(x+h)-f(x-h)}{2h}$	$O(h^2)$

再在$x_{i-2}$、$x_{i+2}$处展开：

$$f_{i+2} = f_i + 2h \cdot f_i' + \frac{4h^2}{2!} f_i'' + \frac{8h^3}{3!} f_i'''+...$$

$$f_{i-2} = f_i - 2h \cdot f_i' + \frac{4h^2}{2!} f_i'' - \frac{8h^3}{3!} f_i'''+...$$

又 $\because$

$$f_{i+1} = f_i + h f'_i + \frac{h^2}{2} f''_i + \frac{h^3}{6} f'''_i + O(h^4)$$

$$f_{i-1} = f_i - h f'_i + \frac{h^2}{2} f''_i - \frac{h^3}{6} f'''_i + O(h^4)$$

消去$f_i'''$，可得：

数值积分

数据滤波