计算物理复习[ep.1|方程数值解法、函数的插值与拟合]

本文含有大量数学公式，请使用电脑或平板电脑浏览以获得最佳阅读体验

Intro

一年两度的期末考试又要来啦~

这篇博客是我大三下学期计算物理课程的复习整理，共包含十个章节，代码部分均以 Matlab 语言写就。参考教材是刘金远等人编著的《计算物理学》，科学出版社的。

复习和平时作业所写的相关代码已经上传到我的 GitHub 仓库中了，仅供学习交流使用。

第一章 方程数值解法

线性方程组

高斯消元法

高斯消元法属于直接法，即：通过有限步运算，减少或消去未知量个数，求得方程精确解，其核心思想为把方程组的系数矩阵化为上三角型，从最后开始，依次求 $x_n,\cdots,x_2,x_1$ 的数值：

$$\begin{array}{c} \left\{\begin{array}{c} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1 \\ a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2 \\ \cdots \\ a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n=b_n \end{array}\right . \\ \Downarrow \\ \left [ \begin{array}{c|c} \begin{matrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \vdots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \\ \end{matrix}& \begin{matrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{matrix} \end{array} \right ] \\ \Downarrow \\ \left [ \begin{array}{c|c} \begin{matrix} a_{11}^0 & a_{12}^0 & \cdots & a_{1n}^0 \\ 0 & a_{22}^1 & \cdots & a_{2n}^1 \\ 0 & 0 & \ddots & \vdots \\ 0 & 0 & \cdots & a_{nn}^{n-1} \\ \end{matrix}& \begin{matrix} b_1^0 \\ b_2^1 \\ \vdots \\ b_n^{n-1} \end{matrix} \end{array} \right ] \end{array}$$

最下面的矩阵中元素的上标表示其被改变的次数

追赶法

追赶法是高斯消去法用于解三对角矩阵的简化形式（因此也属于直接法），这种方法比较简单，计算量小，节省存贮量。

该方法之所以能减少存储量是因为对于一个 n 阶的系数矩阵来说，我们只需要四个长度为 n 的向量来分别存储其三个对角向量 $a_n, b_n, c_n$ 和最右边的 $d_n$：

虽然在矩阵中没有 $a_1$ 和 $c_n$ 元素，但是为了形式上保持一致，我们可以令它们等于零来将 a，c 向量凑够 n 项

$$\begin{bmatrix} b_1 & c_1 & 0 &\cdots & 0 \\ a_2 & b_2 & c_2 &\cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & a_{n-1} & b_{n-1} &c_{n-1} \\ 0 & \cdots & 0 & a_n & b_{n} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_{n-1} \\ x_n \end{bmatrix} = \begin{bmatrix} d_1 \\ d_2 \\ \vdots \\ d_{n-1} \\ d_n \end{bmatrix}$$

和高斯法相似，追赶法将通过一系列运算将 a 向量化为零，b 向量化为1来构造一个只有两个对角向量的 “上三角矩阵”：

$$\left [ \begin{array}{c|c} \begin{matrix} 1 & w_1 & 0 &\cdots & 0 \\ 0 & 1 & w_2 &\cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & 1 &w_{n-1} \\ 0 & \cdots & 0 & 0 & 1 \end{matrix}& \begin{matrix} g_1 \\ g_2 \\ \vdots \\ g_{n-1} \\ g_n \end{matrix} \end{array} \right ]$$

其元素的递推公式如下：

$$\begin{cases} w_1=c_1/b_1\\ w_i=c_i/(b_i-a_iw_{i-1})&(i=2,3,\cdots,n)\\ g_1=d_1/b_1\\ g_i=(d_i-a_ig_{i-1})/(b_i-a_iw_{i-1})&(i=2,3,\cdots,n) \end{cases} \tag{1}$$

最终结果：

$$\begin{cases}x_n=g_n\\ x_i=g_i-w_i x_{i+1}\left(i=n-1,n-2,\cdots1\right)\end{cases} \tag{2}$$

下面是一个简单的示例程序：

%追赶法解方程组
clear

a = [0 -1 -1 -1 -1 -1];
b = [2 2 2 2 2 2];
c = [-1 -1 -1 -1 -1 0];
d = [1 0 1 0 0 1];

n = 6;
w = zeros(1, n);
g = zeros(1, n);

w(1) = c(1)/b(1);
g(1) = d(1)/b(1);

for i=2:n
    w(i) = c(i)/(b(i)-a(i)*w(i-1));
    g(i) = (d(i)-a(i)*g(i-1))/(b(i)-a(i)*w(i-1));
end

x = zeros(1,n);
x(n) = g(n);
for i = n-1:-1:1
    x(i) = g(i)-w(i)*x(i+1);
end

和高斯法一样，追赶法也着有直接法的缺陷：

• 破坏了矩阵本身的系数结构：大型稀疏矩阵，对称阵
• 时间复杂度高，内存消耗大

而迭代法可以在不改变矩阵系数的情况下完成求解。

雅可比迭代法

雅可比法（Jacobi Method）是一种解对角元素几乎都是各行和各列的绝对值最大的值的线性方程组的算法。

保证收敛条件：

系数矩阵为严格或不可约地对角占优矩阵 $|a_{ii}|>\sum\limits_{j\neq i}|a_{ij}|$ （其实就是对角元素是各行和各列的绝对值最大的值）

构造雅可比迭代法的迭代公式步骤如下：

1. 分离对角元

   对于系数矩阵 $A$，有 $A=D+R$：

   $$D=\begin{bmatrix}a_{11}&0&\cdots&0\\ 0&a_{22}&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&a_{nn}\end{bmatrix},\ \ R=\begin{bmatrix}0&a_{12}&\cdots&a_{1n}\\ a_{21}&0&\cdots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{n1}&a_{n2}&\cdots&0\end{bmatrix}$$

2. 改写方程组

   $$\mathbf{Ax}=\mathbf{b}\Leftrightarrow(\mathbf{D}+\mathbf{R})\mathbf{x}=\mathbf{b}\Leftrightarrow\mathbf{D}\mathbf{x}=\mathbf{b}-\mathbf{R}\mathbf{x}\Leftrightarrow\mathbf{x}=\mathbf{D}^{-1}(\mathbf{b}-\mathbf{R}\mathbf{x}) \tag{3}$$

3. 构造不动点迭代方程

   $$\mathbf x_{k+1}=\mathbf D^{-1}(\mathbf b-\mathbf R\mathbf x_k) \tag{4}$$

   其中各个元素可以表示为:

   $$x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum_{j\neq i}a_{ij}x_j^{(k)}\right),\quad i=1,2,\ldots,n. \tag{5}$$

下面这个实例可以帮助你理解：

$$\begin{cases}10x_1-x_2-2x_3=7.2\\ -x_1+10x_2-2x_3=8.3\\ -x_1-x_2+5x_3=4.2\end{cases} \ \Rightarrow\ \begin{cases}x_1=0.1x_2+0.2x_3+0.72\\ x_2=0.1x_1+0.2x_3+0.83\\ x_3=0.2x_1+0.2x_2+0.84\end{cases} \tag{6}$$

你会发现，上面那个步骤说了一大堆其实就想表达一个思想：把每行占优系数所对应的变量单独提出来放到等号左边。。。（明明一句话的事情教科书上非要搞这么麻烦）

高斯-赛德尔迭代法

高斯-赛德尔（Gauss-Seidel）迭代法（以下简称为 GS 迭代法），其实就是雅可比迭代法的优化版本（优化不多，但有用）。

以迭代公式（6）来举例说明雅可比迭代法和 GS 迭代法的区别：

雅可比迭代法：

$$\begin{cases} x_1^{k+1}=0.1x_2^{k}+0.2x_3^{k}+0.72\\ x_2^{k+1}=0.1x_1^{k}+0.2x_3^{k}+0.83\\ x_3^{k+1}=0.2x_1^{k}+0.2x_2^{k}+0.84 \end{cases} \tag{7}$$

GS 迭代法：

$$\begin{cases} x_1^{k+1}=0.1x_2^{k}+0.2x_3^{k}+0.72\\ x_2^{k+1}=0.1x_1^{k+1}+0.2x_3^{k}+0.83\\ x_3^{k+1}=0.2x_1^{k+1}+0.2x_2^{k+1}+0.84 \end{cases} \tag{8}$$

其中 $k$ 为迭代次数。怎么样，发现不同之处了吗？雅可比迭代法每次迭代得到的 $\mathbf{x}^{k+1}$ 全部用的都是上一次迭代得到的 $\mathbf{x}^{k}(x_1^{k},x_2^{k},x_3^{k})$ ，而 GS 迭代法在每求得一个 $x_n^{k}$ 之后都会用于求下一个 $x_{n+1}^{k}$ ，换句话说，就是 GS 迭代法”更新“的速度比雅可比要快，这一点优化可以让 GS 迭代法更快地收敛。

非线性方程

二分法

二分法适用于明确某区间 [a,b] 含有根的情况，其运行原理就是通过不断缩小含根区间（每次都砍半，这也是它名字的来源）来确定不同精度下的根的大致范围。对二分法来说，决定是否继续缩小区间的判据为现有区间大小是否满足精度要求，而决定选择哪一边继续缩小则取决于区间两端的函数值是否符号相反。

下面是一个简单的示例代码：

%二分法求方程根
a = 2;
b = 1.5;
eps = 1e-5;
err = 10;

while (err > eps)
    x = (a + b)/2;
    if f(x)*f(a) < 0
        b = x;
    else
        a = x;
    end
    err = abs(f(x));
end

function y = f(x)
    y = sin(x)-x^2/4;
end

不动点迭代收敛判据

• 判据1：设 $x^*=g(x^*)$， 若 $g(x)$ 在 $x^*$ 附近连续可微且 $|g^{\prime}(x^{*})|<1$，则迭代格式 $x_{k+1}=g(x_k)$ 在 $x^*$ 附近局部收敛。

  注: 由于 $x^*$ 事先未知，故实际应用时，代之以近似判则 $\big|g'(x_0)\big|<1$。但需注意，这实际上是假设了 $x_0$ 充分接近 $x^*$ ，若 $x_0$ 离 $x^*$ 较远，迭代格式可能不收敛。

• 判据2 (非局部收敛定理）如果 在 $[a,b]$ 上连续可微且以下条件满足：

  1. 若 $x\in[a,b]$，则 $g(x)\in[a,b]$
  2. 对 $\forall x\in[a,b],|g'(x)|\le L<1$

  那么，对 $\forall x\in[a,b],x_{k+1}=g(x_k)$ 的解序列 $\{x_k\}$ 收敛于 $x^*$

代码判敛

用代码进行判敛就方便多了，只需要设定精度（eps）和初始差值（err），我们就可以将本次迭代结果 y 和上一次的结果 x 的范数（norm）相减得到差值，并将其和精度进行对比来决定是否继续迭代，代码如下：

eps = 1e-5;
err = 10;
x = [0.1 0.1 -0.1];
while err > eps
    y = g(x);
    err = norm(y - x);
    x = y;
end

牛顿迭代法

线性方程容易解，但是它终归只占了世间千万方程的一小小小部分，我们在工作学习中遇到的绝大多数函数或者方程都是非线性的，这就对求解造成了很大的困扰。但是别忘了，我们对付非线性方程也有一个大杀器——泰勒展开。如果函数 $f(x)$ 是连续可微的， $x^*$ 是 $f(x)=0$ 的实根，$x_k$ 是迭代过程中的近似根，那我们就可以利用泰勒展开把 $f(x)$ 在 $x_k$ 附近展开为泰勒级数，然后只取其线性部分进行计算。这也是牛顿迭代法的基本思想。

设 $f(x)$ 在其零点 $x^*$ 附近连续可微，已知 $x_k$ 为的第 k 次近似值，则：

$$\begin{array}{l} f(x^*)=f(x_k+x^*-x_k)=f(x_k)+f'(x_k)(x^*-x_k)+\mathcal{O}\Big((x^*-x_k)^2\Big) \\ f(x^*) =0\to f(x_k)+f'(x_k)(x^*-x_k)\approx0 \end{array} \tag{9}$$

构造不动点迭代式：

$$\begin{gathered} f(x_k)+f'(x_k)(x_{k+1}-x_k)=0 \\ x_{k+1}=x_{k}-{\frac{f(x_{k})}{f^{\prime}(x_{k})}} \end{gathered} \tag{10}$$

得到迭代函数：

$$g(x)=x-{\frac{f(x)}{f^{\prime}(x)}} \tag{11}$$

非线性方程组

• 非线性方程组的数值解通常采用迭代法
• 雅可比迭代及赛德尔迭代，均可用于求解非线性方程组，求解单一非线性方程的牛顿法和弦截法等，也可推广到求解非线性方程组。

雅可比和 GS 迭代法用于解非线性方程组

以求解下面非线性方程组为例：

$$\begin{cases}3x_1-\cos(x_2x_3)-\frac{1}{2}=0\\ x_1^2-81(x_2+0.1)^2+\sin x_3+1.06=0\\ \exp(-x_1x_2)+20x_3+\frac{(10\pi-3)}{3}=0\end{cases} \tag{12}$$

取初值 $\mathbf{x^0}=(0.1,0.1,-0.1)^T$，精度（eps）取 $10^{-5}$

虽然换成了非线性方程组，但是基本思想没变，即：每一行分别提出来一个 $x_n$ ，而且行与行的 $x_n$ 互不相同（因为是非线性方程组，所以此处的”占优项“概念变更为：容易提出 $x_n$ 的项）

下式是一种构造结果：

$$\left \{\begin{aligned} &x_1 =\frac{1}{3}\cos(x_2x_3)+\frac{1}{6} \\ &x_2 =\frac{1}{9}(x_1^2+\sin x_3+1.06)^{0.5}-0.1 \\ &x_3 =-\exp(-x_1x_2)/20-\frac{(10\pi-3)}{60} \end{aligned}\right . \tag{13}$$

雅可比迭代法方程组：

function y = g(x)
    y(1) = 1/3*cos(x(2)*x(3)) + 1/6;
    y(2) = 1/9*(x(1)^2 + sin(x(3)) + 1.06)^0.5 - 0.1;
    y(3) = -1*exp(-x(1)*x(2))/20 - (10*pi - 3)/60;
end

GS 迭代法方程组：

function y = g(x)
    y(1) = 1/3*cos(x(2)*x(3)) + 1/6;
    y(2) = 1/9*(y(1)^2 + sin(x(3)) + 1.06)^0.5 - 0.1;
    y(3) = -1*exp(-y(1)*y(2))/20 - (10*pi - 3)/60;
end

从代码也可以看出雅可比迭代法和 GS 迭代法的区别，即 GS 迭代法在求 y(2) 时用到了 y(1) 并非雅可比迭代法中的 x(1) ，求 y(3) 时也是如此。

牛顿迭代法解非线性方程组

对于非线性方程组 $F(x) = 0$，牛顿迭代法的形式为：

$$\mathrm{x_{k+1}=x_k-J^{-1}(x_k)F(x_k)} \tag{14}$$

其中 $J$ 为非线性方程组对应的 Jacobi 矩阵（非线性方程求导数，非线性方程组求 Jacobi 矩阵），下面将以此为根据构造迭代函数。

由 Jacobi 矩阵的定义：

$$J = \begin{bmatrix}\dfrac{\partial f_1}{\partial x_1}&\cdots&\dfrac{\partial f_1}{\partial x_n}\\ \vdots&\ddots&\vdots\\ \dfrac{\partial f_m}{\partial x_1}&\cdots&\dfrac{\partial f_m}{\partial x_n}\end{bmatrix} \tag{15}$$

其中矩阵元可表示为：

$$J_{ij} = \frac{\partial f_i}{\partial x_j} \tag{16}$$

对于本题：

$$F(x) = \begin{bmatrix} f_1\\ f_2\\ f_3 \end{bmatrix} =\begin{bmatrix} 3x_1 - \mathrm{cos}(x_2x_3)-\frac{1}{2} \\ x_1^2-81(x_2+0.1)^2+\mathrm{sin}x_3+1.06 \\ e^{-x_1x_2}+20x_3+\frac{10\pi - 3}{3} \end{bmatrix} \tag{17}$$

将（17）代入（15）得到和 $F(x)$ 对应的 Jacobi 矩阵：

$$J = \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \frac{\partial f_1}{\partial x_3}\\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \frac{\partial f_2}{\partial x_3}\\ \frac{\partial f_3}{\partial x_1} & \frac{\partial f_3}{\partial x_2} & \frac{\partial f_3}{\partial x_3}\\ \end{bmatrix} = \begin{bmatrix} 3 & x_3\mathrm{sin}(x_2x_3) & x_2\mathrm{sin}(x_2x_3) \\ 2x_1 & -162(x_2+0.1) & \mathrm{cos}x_3 \\ -x_2e^{-x_1x_2} & -x_1e^{-x_1x_2} & 20 \end{bmatrix} \tag{18}$$

抱歉我的 LaTeX 渲染器太拉了，把这个👆 Jacobi 矩阵搞得好难看

由式（14）得到迭代函数：

$$g(x) = x - J^{-1}F$$

其中的 $J^{-1}F$ 用 Matlab 的 J\F 计算来代替：

function y = g(x)           %迭代函数
    y = x - J(x)\F(x);
end

第二章 函数的插值与拟合

插值和拟合的区别

如图，左侧为插值结果，右侧为拟合结果

插值和拟合的最大不同之处为：插值要求目标函数必须经过给定的数据点，而拟合只表达数据的变化趋势；插值和拟合的另一个不同点为：插值方法反映了给定一组分立数据的局域性质（分段插值），而拟合方法给出了这组数据的整体性质（全部数据）

对于一组给定数据 $\mathcal{D}=\{(x_1,y_1),(x_2,y_2),\ldots,(x_k,y_k)\}$ ，如果我们可以构造一个简单函数 $\varphi(x)$ 作为函数 $y = f(x)$ 的近似表达式，即 $y=f(x)\approx\varphi(x)$ ，使得：

$$\varphi(x_1)=y_1,\varphi(x_2)=y_2,...,\varphi(x_k)=y_k$$

这类问题称为插值问题。 $f(x)$ 称为被插值函数，$\varphi(x)$ 称为插值函数， $x_1 ,x_2, ... , x_k$ 称为插值节点。

多项式插值

一种简单的插值方法，思路是使用多项式函数：

$$P_n(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n \tag{19}$$

来逼近被插值函数 $f(x)$ ，我们可以列出一组 k 个方程式：

$$\begin{cases} P_n(x_1)=a_0+a_1x_1+a_2x_1^2+\cdots+a_nx_1^n=y_1 \\ P_n(x_2)=a_0+a_1x_2+a_2x_2^2+\cdots+a_nx_2^n=y_2 \\ \vdots \\ P_n(x_k)=a_0+a_1x_k+a_2x_k^2+\cdots+a_nx_k^n=y_k \end{cases} \tag{20}$$

由于是多项式，式（20）的系数矩阵为范德蒙行列式：

$$\begin{bmatrix}1&x_1&x_1^2&\cdots&x_1^n\\ 1&x_2&x_2^2&\cdots&x_2^n\\ 1&x_3&x_3^2&\cdots&x_3^n\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&x_k&x_k^2&\cdots&x_k^n\end{bmatrix}\begin{bmatrix}a_0\\ a_1\\ a_2\\ \vdots\\ a_n\end{bmatrix}=\begin{bmatrix}y_1\\ y_2\\ \vdots\\ y_3\\ y_k\end{bmatrix} \tag{21}$$

我所剩不多的线性代数知识告诉我，未知数的数量要大于等于方程的数量才能保证有解（不一定是唯一解），即 $n+1 \ge k$ ，因此多项式函数的阶数（n）必须大于或等于数据点个数（k）-1

拉格朗日插值法

让我们先看看对拉格朗日基函数数学上的定义：

对给定数据 $\mathcal{D}=\{(x_1,y_1),(x_2,y_2),\ldots,(x_k,y_k)\}$ ，有插值函数

$$P_n(x)=y_1l_1(x)+y_2l_2(x)+\cdots+y_kl_k(x) \tag{22}$$

其中的基函数 $l_i(x)$ 满足：

$$\left\{\begin{array}{l}{l_{i}(x_{i})=1}\\ {l_{i}(x_{j})=0(j\neq i)}\\ \end{array}\right. \ \Rightarrow\ \begin{cases} P_n(x_1)=y_1\\ P_n(x_2) =y_{2}\\ \vdots\\ P_n(x_k)=y_k \end{cases} \tag{23}$$

让我们先考虑最简单的情况，只有两个点 $\mathcal{D}=\{(x_1,y_1),(x_2,y_2)\}$ ，这种情况对应的插值函数是什么样的呢？很明显，在各种各样的函数中，一元一次函数是既符合条件又很简单的理想函数，毕竟”两点可确定一条直线“，于是，我们就有了：

$$\left\{\begin{array}{l} {l_{1}(x_{1})=1}\\ {l_{1}(x_{2})=0}\\ \end{array}\right. \ \Rightarrow\ l_1(x)=\frac{x-x_2}{x_1-x_2} \tag{24}$$

同理：

$$\left\{\begin{array}{l} {l_{2}(x_{2})=1}\\ {l_{2}(x_{1})=0}\\ \end{array}\right. \ \Rightarrow\ l_2(x)=\frac{x-x_1}{x_2-x_1} \tag{24}$$

如果是三个点 $\mathcal{D}=\{(x_1,y_1),(x_2,y_2),(x_3,y_3)\}$ 呢？对第一个基函数 $l_1(x)$ 来说我们需要构建一个同时满足：

$$\begin{matrix} \begin{cases} {l_{1}(x_{1})=1}\\ {l_{1}(x_{2})=0}\\ \end{cases} & & \begin{cases} {l_{1}(x_{1})=1}\\ {l_{1}(x_{3})=0}\\ \end{cases} \end{matrix}$$

很明显，这两个条件是”与“的关系，因此我们可以用乘法将它们联系起来：

$$\begin{array}{c} \begin{matrix} \begin{cases} {l_{1}(x_{1})=1}\\ {l_{1}(x_{2})=0}\\ \end{cases} & & \begin{cases} {l_{1}(x_{1})=1}\\ {l_{1}(x_{3})=0}\\ \end{cases} \\ \Downarrow & & \Downarrow \\ \frac{x-x_2}{x_1-x_2} & \times & \frac{x-x_3}{x_1-x_3} \\ & \Downarrow & \end{matrix} \\ l_1(x)=\frac{x-x_2}{x_1-x_2} \cdot \frac{x-x_3}{x_1-x_3} \end{array} \tag{25}$$

以此类推，对 $\mathcal{D}=\{(x_1,y_1),(x_2,y_2),\ldots,(x_k,y_k)\}$ ，第 $i$ 个拉格朗日基函数为：

$$\begin{array}{c} l_i(x)=\prod_{p=2}^k\frac{x-x_p}{x_i-x_p},i=1\\ l_i(x)=\prod_{p=1,p \ne i}^k\frac{x-x_p}{x_i-x_p},i\ne1 \end{array} \tag{26}$$

牛顿插值法&差商表

拉格朗日插值虽好可不要贪杯哦~，但是由于 $l(x)$ 依赖于全部基点，若计算出所有基函数之后又要增加数据点，则全部 $l(x)$ 都需要重算。为避免以上麻烦，我们可以使用牛顿插值法。

牛顿插值基函数定义：$\phi_{0}(x)=1,\phi_{i}(x)=(x-x_{i})\phi_{i-1}(x)\quad(1\leq i\leq k)$

从上面的定义来看，我们会发现，在增加数据点 $x_i$ 时，牛顿插值法的基函数可以通过在已有的基函数之上乘上一个 $(x-x_i)$ 得到，无需全部重算，此时插值函数可以表示为：

$$N_k(x)=c_0\phi_0+c_1\phi_1+\cdots+c_k\phi_k \tag{27}$$

通过 $\mathcal{D}=\{(x_1,y_1),(x_2,y_2),\ldots,(x_k,y_k)\}$ 我们可以求各项系数 $c_i$

$$\begin{array}{c} \begin{array}{l} N_{k}(x_{1}) =y_1=c_0 \\ N_{k}(x_{2}) =y_2=c_0+c_1(x_2-x_1) \\ N_{k}(x_{3}) =y_3=c_0+c_1(x_3-x_1)+c_2(x_3-x_1)(x_3-x_2) \\ \vdots \\ N_{k}(x_{k})=y_{k}=c_{0}+c_{1}(x_{k}-x_{1})+\cdots+c_{k-1}(x_{k}-x_{1})(x_{n}-x_{2})\cdots(x_{k}-x_{k-1}) \end{array}\\ \\ \Downarrow \\ \begin{aligned} &c_{0} =y_1 \\ &c_{1} =\frac{y_2-y_1}{x_2-x_1} \\ &c_{2}=\left({\frac{y_{3}-y_{2}}{x_{3}-x_{2}}}-{\frac{y_{2}-y_{1}}{x_{2}-x_{1}}}\right)/(x_{3}-x_{1}) \\ \vdots \end{aligned} \end{array}$$

对于一个数据表 $\{(x_1,f(x_1)),(x_2,f(x_2)),\ldots,(x_k,f(x_k))\}$

我们将 $f[x_i]$ 定义为 $f(x)$ 关于 $x_i$ 的零阶差商， $f(x)$ 关于 $x_i,x_{i+1}$ 的一阶差商定义为：

$$f[x_i,x_{i+1}]=\frac{f[x_{i+1}]-f[x_i]}{x_{i+1}-x_i}$$

推而广之， $f(x)$ 关于 $x_i,x_{i+1},\cdots,x_{i+m}$ 的 $m$ 阶差商定义为：

$$f[x_i,x_{i+1},\ldots,x_{i+m}]=\frac{f[x_{i+1},x_{i+2},\ldots,x_{i+m}]-f[x_i,x_{i+1},\ldots,x_{i+m-1}]}{x_{i+m}-x_i} \tag{28}$$

和牛顿插值法的第 $i$ 项对比：

$$n_i(x)=f[x_1,x_2,\dots,x_{i+1}]\phi_i(x) \tag{29}$$

我们得知，牛顿插值法的插值函数系数 $c_i$ 为 $f(x)$ 关于 $x_1,x_{2},\cdots,x_{i+1}$ 的 $i$ 阶差商

作下图差商表：

上表对角项 $y_1,f[x_1,x_2],\cdots,f[x_1,x_2,\ldots,x_{k}]$ 就对应着 $c_1,c_2,\cdots,c_{k}$

注意，牛顿插值法的插值函数 $N_k$ 中的 $k$ 和数据点要根据实际情况选取，如果你要预测的点在 $x_3$ 和 $x_4$ 之间而且要求 $k=2$ ，那么利用 $\mathcal{D}=\{(x_2,y_2),(x_3,y_3),(x_4,y_4)\}$ 或 $\mathcal{D}=\{(x_3,y_3),(x_4,y_4),(x_5,y_5)\}$ 进行插值是不错的选择

三次样条插值

所谓的三次样条插值算法，让我用下图中这些波浪线和点为例来解释 我知道可能有点简陋，您嘞将就着看吧

${\Huge \cdot\sim\cdot\sim\cdot\sim\cdot\sim{\normalsize\cdots}\sim\cdot\sim\cdot\sim\cdot\sim\cdot\sim\cdot}$

这些点代表着给定的数据集 $\mathcal{D}=\{(x_1,y_1),(x_2,y_2),\ldots,(x_k,y_k)\}$ ，而中间的波浪线就代表着每一段插值函数 $S_i(x)$，总的插值函数 $S(x)$ 应该满足：

$$S(x)=\left\{\begin{matrix}S_1(x),&x\in[x_1,x_2]\\ S_2(x),&x\in[x_2,x_3]\\ \vdots&\vdots\\ S_{k-1}(x),&x\in[x_{k-1},x_k]\end{matrix}\right. \space \space \space \space \space \space \space \space \space \space S_i(x)=\frac{1}{6}a_ix^3+\frac{1}{2}b_ix^2+c_ix+d_i \tag{30}$$

很明显，$S(x)$ 是一个分段函数，而且其每一段都是一个三次函数（毕竟是"三次"样条插值嘛），同时也别忘了，$S(x)$ 作为一个插值函数，要满足插值条件和连续条件，i.e. $S(x)$ 经过所有给定的数据点并且在定义域 $[x_1,x_k]$ 内具有连续的一阶和二阶导数，用数学公式来解释就是：

$$\left\{\begin{array}{l l} {S_{i-1}(x_{i})=S_{i}(x_{i})}&{i=2,3,\cdots,k-1}\\ {S_{i-1}^{\prime}(x_{i})=S_{i}^{\prime}(x_{i})}&{i=2,3,\cdots,k-1}\\ {S_{i-1}^{\prime\prime}(x_{i})=S_{i}^{\prime\prime}(x_{i})}&{i=2,3,\cdots,k-1} \\ S(x_1)=y_1,S(x_2)=y_2,\cdots,S(x_k)=y_k \end{array}\right. \tag{31}$$

推公式预警！看不下去的可以直接去看小结😀

好了，我们现在需要求解 $k-1$ 个分段的插值函数 $S_i(x)$ ，它们每一个都有四个待定的系数，因此总共要解出 $4(k-1)=4k-4$ 个系数。但是根据（31），一共只有 $3(k-2)+k=4k-6$ 个约束条件，还差了两个，我们就要通过引入边界条件来解决这个问题

现在我们求 $S_i(x)$ 函数，先介绍补充两端的二阶导数值（三弯矩算法）。我们无法直接得到三次函数 $S_i(x)$ ，但是我们可以先挑个软的柿子捏，即 $S_i(x)$ 的二阶导数 $S_i''(x)$ ，因为是一次函数所以很适合作为求解的起点，令：

$$S_j''(x_j)=M_j\space(j=1,2,3,\cdots,k),\space h_j=x_{j+1}-x_j \tag{32}$$

因为 $S_j''(x_j)$ 过 $(x_j,M_j),(x_j,M_{j+1})$ ，有：

$$S_j''(x)=\frac{x_{j+1}-x}{h_j}M_j+\frac{x-x_j}{h_j}M_{j+1} \tag{33}$$

积两次分：

$$S_{j}(x)={\frac{(x_{j+1}-x)^{3}}{6h_{j}}}M_{j}+{\frac{(x-x_{j})^{3}}{6h_{j}}}M_{j+1}+c_{1}x+c_{2} \tag{34}$$

联立 $\begin{cases}S_j(x_j)=\frac{1}{6}h_j^2M_j+c_1x_j+c_2=y_j\\ S_j(x_{j+1})=\frac{1}{6}h_j^2M_{j+1}+c_1x_{j+1}+c_2=y_{j+1}\end{cases}$ 解出来 $c_1$ 和 $c_2$ 再回代入（34）得：

$$\begin{array}{c} S_j(x)=\frac{(x_{j+1}-x)^3}{6h_j}M_j+\frac{(x-x_j)^3}{6h_j}M_{j+1}+\left(y_j-\frac{M_jh_j^2}{6}\right)\frac{x_{j+1}-x}{h_j}+\left(y_{j+1}-\frac{M_{j+1}h_j^2}{6}\right)\frac{x-x_j}{h_j} \\ \Downarrow \space \mathrm{Derivative} \\ S_j'(x)=-\frac{(x_{j+1}-x)^2}{2h_j}M_j+\frac{(x-x_j)^2}{2h_j}M_{j+1}+\frac{y_{j+1}-y_j}{h_j}-\frac{M_{j+1}-M_j}{6}h_j \end{array} \tag{35}$$

由于 $S_{j}^{\prime}(x_{j})=S_{j-1}^{\prime}(x_{j})$ 有：

$$-\frac{h_j}{3}M_j-\frac{h_j}{6}M_{j+1}+\frac{y_{j+1}-y_j}{h_j}=\frac{h_{j-1}}{6}M_{j-1}+\frac{h_{j-1}}{3}M_j+\frac{y_j-y_{j-1}}{h_{j-1}}$$

化简得到：

$${\frac{h_{j-1}}{h_{j-1}+h_{j}}}M_{j-1}+2M_{j}+{\frac{h_{j}}{h_{j-1}+h_{j}}}M_{j+1}={\frac{6}{h_{j-1}+h_{j}}}\left({\frac{y_{j+1}-y_{j}}{h_{j}}}-{\frac{y_{j}-y_{j-1}}{h_{j-1}}}\right) \tag{36}$$

令 $\mu_{j}={\frac{h_{j-1}}{h_{j-1}+h_{j}}},\lambda_{j}={\frac{h_{j}}{h_{j-1}+h_{j}}},d_{j}={\frac{6}{h_{j-1}+h_{j}}}\Biggl[{\frac{y_{j+1}-y_{j}}{h_{j}}}-{\frac{y_{j}-y_{j-1}}{h_{j-1}}}\Biggr]$ 有：

$$\begin{array}{c} \mu_j M_{j-1}+2M_j+\lambda_j M_{j+1}=d_j\quad j=2,3,\cdots,k-1 \\ \Downarrow \\ \begin{cases}{\mu_{2}M_{1}+2M_{2}+\lambda_{2}M_{3}=d_{2}}\\ {\mu_{3}M_{2}+2M_{3}+\lambda_{3}M_{4}=d_{3}}\\ {\ldots}\\ {\mu_{k-2}M_{k-3}+2M_{k-2}+\lambda_{k-2}M_{k-1}=d_{k-2}}\\ {\mu_{k-1}M_{k-2}+2M_{k-1}+\lambda_{k-1}M_{k}=d_{k-1}}\\ \end{cases} \end{array} \tag{37}$$

注意，$M_1$ 和 $M_k$ 是作为边界条件补充进来的，因此它们是已知的，做适当移项并用矩阵表示上方程组得：

$$\begin{bmatrix} 2&\lambda_2\\ \mu_3&2&\lambda_3\\ &\ddots&\ddots&\ddots\\ &&\mu_{k-2}&2&\lambda_{k-2}\\ &&&\mu_{k-1}&2 \end{bmatrix} \begin{bmatrix} M_2\\ M_3\\ \vdots\\ M_{k-2}\\M_{k-1} \end{bmatrix}= \begin{bmatrix} d_2-\mu_2M_1\\ d_3\\ \vdots\\ d_{k-2}\\ d_{k-1}-\lambda_{k-1}M_k \end{bmatrix}$$

至此，我们就把求 $k-1$ 个分段的插值函数 $S_i(x)$ 的问题变成了求解严格对角占优三对角矩阵问题，用追赶法或者其他迭代法求解即可。

有人会问，你讲了这么大一堆都是在用二阶导数的边界条件，那一阶的怎么办捏？

莫急，从二阶换到一阶其实很简单，由（35）中 $S_j'(x)$ 的表达式和边界条件 $\left\{\begin{array}{l l}{S^{\prime}(x_{1})=S_{1}^{\prime}(x_{1})=m_{1}}\\ {S^{\prime}(x_{k})=S_{k-1}^{\prime}(x_{k})=m_{k}}\\ \end{array}\right.$ 可以得到：

$$\begin{aligned} &S_{1}^{\prime}(x_{1})=-\frac{h_{1}}{3}M_{1}-\frac{h_{1}}{6}M_{2}+\frac{y_{2}-y_{1}}{h_{1}}=m_{1} \\ &S_{k-1}^{\prime}(x_{k})=\frac{h_{k-1}}{6}M_{k-1}+\frac{h_{k-1}}{3}M_{k}+\frac{y_{k}-y_{k-1}}{h_{k-1}}=m_{k} \end{aligned} \tag{38}$$

把 $M_1,M_k$ 反解出来：

$$\begin{aligned} &M_1=-\frac{3}{h_1}m_1-\frac{1}{2}M_2+\frac{3}{h_1^2}(y_2-y_1) \\ &M_k=\frac{3}{h_{k-1}}m_k-\frac{1}{2}M_{k-1}-\frac{3}{h_{k-1}^2}(y_k-y_{k-1}) \end{aligned} \tag{39}$$

再代入（37）即可得到和一阶边界条件对应的方程组

三次样条插值小结

对于一组数据 $\mathcal{D}=\{(x_1,y_1),(x_2,y_2),\ldots,(x_k,y_k)\}$

如果用三弯矩法，第二边界条件给定：

$$\begin{cases}S''(x_1)=S''_1(x_1)=M_1\\ S''(x_k)=S''_{k-1}(x_k)=M_k\end{cases}$$

有：$S_j(x)=\frac{(x_{j+1}-x)^3}{6h_j}M_j+\frac{(x-x_j)^3}{6h_j}M_{j+1}+\left(y_j-\frac{M_jh_j^2}{6}\right)\frac{x_{j+1}-x}{h_j}+\left(y_{j+1}-\frac{M_{j+1}h_j^2}{6}\right)\frac{x-x_j}{h_j}$

解系数矩阵求 $M_j(j\ne1,k)$：

$$\begin{bmatrix} 2&\lambda_2\\ \mu_3&2&\lambda_3\\ &\ddots&\ddots&\ddots\\ &&\mu_{k-2}&2&\lambda_{k-2} \end{bmatrix} \begin{bmatrix} M_2\\ M_3\\ \vdots\\ M_{k-2} \end{bmatrix}= \begin{bmatrix} d_2-\mu_2M_1\\ d_3\\ \vdots\\ d_{k-2}\\ d_{k-1}-\lambda_{k-1}M_3 \end{bmatrix}$$

其中：

$$\begin{cases} h_j=x_{j+1}-x_j \\ \mu_{j}={\frac{h_{j-1}}{h_{j-1}+h_{j}}}\\ \lambda_{j}={\frac{h_{j}}{h_{j-1}+h_{j}}}\\ d_{j}={\frac{6}{h_{j-1}+h_{j}}}({\frac{y_{j+1}-y_{j}}{h_{j}}}-{\frac{y_{j}-y_{j-1}}{h_{j-1}}}) \end{cases}$$

用第一边界条件 $\left\{\begin{array}{l l}{S^{\prime}(x_{1})=S_{1}^{\prime}(x_{1})=m_{1}}\\ {S^{\prime}(x_{k})=S_{k-1}^{\prime}(x_{k})=m_{k}}\\ \end{array}\right.$ 有：

$$\begin{cases} {(2-\frac{\mu_{2}}{2})M_{2}+\lambda_{2}M_{3}=d_{2}-\mu_2(\frac{3}{h_1^2}(y_2-y_1)-\frac{3}{h_1}m_1)}\\ {\mu_{3}M_{2}+2M_{3}+\lambda_{3}M_{4}=d_{3}}\\ {\ldots}\\ {\mu_{k-2}M_{k-3}+2M_{k-2}+\lambda_{k-2}M_{k-1}=d_{k-2}}\\ {\mu_{k-1}M_{k-2}+(2-\frac{\lambda_{k-1}}{2})M_{k-1}=d_{k-1}}-\lambda_{k-1}(\frac{3}{h_{k-1}^2}(y_2-y_1)-\frac{3}{h_{k-1}}m_k)\\ \end{cases}$$

用矩阵表示：

$$\begin{bmatrix} 2-\frac{\mu_{2}}{2}&\lambda_2\\ \mu_3&2&\lambda_3\\ &\ddots&\ddots&\ddots\\ &&\mu_{k-2}&2-\frac{\lambda_{k-1}}{2}&\lambda_{k-2} \end{bmatrix} \begin{bmatrix} M_2\\ M_3\\ \vdots\\ M_{k-2} \end{bmatrix}= \begin{bmatrix} d_{2}-\mu_2(\frac{3}{h_1^2}(y_2-y_1)-\frac{3}{h_1}m_1)\\ d_3\\ \vdots\\ d_{k-2}\\ d_{k-1}-\lambda_{k-1}(\frac{3}{h_{k-1}^2}(y_2-y_1)-\frac{3}{h_{k-1}}m_k) \end{bmatrix}$$

梯度下降算法

其实拟合相关的最重要的算法是回归算法，但是由于这次期末考试不考（真的是因为这个，绝对不是什么偷懒之类的原因哦）俺就没写。

先上图和代码，下面你看到的这张图是利用梯度下降法找 $f(x)=x^2$ 这个函数最低点迭代20次之后的结果

介是代码👇：

% 梯度下降法找y=x^2最低点
f = @(x)(x.^2);
t = -4:0.01:12;
plot(t,f(t));
times = 20;
alpha = 0.1;
x = 10;
hold on
for i = 1:times
    x1 = x;
    y1 = f(x1);
    fprintf('第%d次迭代:x1=%d,y1=%d\n',i,x1,y1)
    scatter(x1,y1,'red')
    x = x - alpha.*2.*x;
end
grid on
title("Gradient Descent")

通过阅读代码我们可以大致了解梯度下降算法到底在做什么，我们为迭代对象设了个初值 x = 10 ，说明我们迭代的起点是10。如果令迭代递推式为 $x_{i+1}=x_i+\delta_i$ ，那么每次迭代后 $x_i$ 的变化量就是我们下一步要关心的东西了，因为它决定了下一次迭代的“方向”和“步长”。为了得到这个变化量 $\delta_i$，我们就要以梯度为核心展开讨论了。

函数上的每一个点都有其相对应的梯度（导数），对于函数 $f(x)=x^2$ 来说，在所求的 (0,0) 点之后的梯度是正数，反之为负数。为了每一次迭代都可以让 $x_i$ 更接近要求的值，我们可以这样定义 $\delta_i$ ：

$$\delta_i(x)=-\alpha f'(x) \tag{40}$$

因为有一个负号在前，所以 $\delta_i$ 和函数的导数 $f'(x)$ 的符号是相反的，这意味着：如果 $x_i$ 在目标值右边，那么它的导数将会使它有向左边“走”的趋势；如果 $x_i$ 在目标值左边，那么它的导数将会使它有向右边“走”的趋势，在不断的迭代过程中，$x_i$ 将始终受到某种“神秘力量”让它只会向目标值的方向走去。

更妙的是，$x_i$ 离我们的目标点越远，其梯度的绝对值 $|f'(x)|=|2x|$ 就会越大，$x_i$ 就会“走”地越快，就像我妈在离家很远的地方突然想起来家里锅中还煮着东西时的样子；而离我们的目标点越近，$x_i$ 就会“走”地越慢，就像你和你的朋友见面时开着手机定位互相找一样，生怕错过对方。（40）式里还有个 $\alpha$ ，这个就是我帮助 $x_i$ 控制它移动速度的参数，本案例中 alpha = 0.1 ，有了它，$x_i$ 就更不容易“跨过”我们要寻找的目标点了。

由于有以上两个限制条件，在一次次迭代过程中，我们的 $x_i$ 会越来越逼近目标点，本案例中是通过 for 循环来控制迭代次数的，其实也可以模仿 2.2.3 节中利用 while 来进行代码判敛，具体实现方式不再赘述。

終了

OK 这篇文章写到这里我觉得足够长了，（向屏幕一角看了一眼）好家伙，Typora 统计一共有8642个字，顶上一篇小论文了（我写了快三天啊，三天啊三天！），感谢你能保持耐心看到这里，为你和我自己都点个赞~

这是计算物理复习的第一篇文章，接下来我计划用总共四篇文章来总结下这门课我认为需要记住的知识点。最近肝一点，争取在三天内全部搞出来吧 Orz

最后，再次感谢你能看完这篇文章，有误之处欢迎反馈！