第 14 章 统计计算方法

CHAPTER STATUS编写进度10%查看全书进度

14.1 非线性方程的求解

14.1.0.1 二分法

Method 14.1. (Bisection Method)

设函数\(f\)\([a,b]\)上连续,且满足\(f(a)f(b)<0\)。为求解非线性方程\(f(x)=0\),由引理 3.3可知,存在至少一个根\(x^{\star}\in(a,b)\)。二分法通过反复对区间进行二等分来逐步缩小根所在的区间:令:

\[\begin{equation*} c=\frac{a+b}{2} \end{equation*}\]

\(f(c)=0\),则\(c\)即为方程的解;否则根据符号判定准则,取满足:

\[\begin{equation*} f(a)f(c)<0\quad\text{或}\quad f(c)f(b)<0 \end{equation*}\]

的子区间作为新的区间\([a,b]\),并重复上述过程。

经过\(n\)次迭代后,所得区间长度不超过\((b-a)/2^n\),任意迭代点\(x_n\)与真根\(x^{\star}\)的误差满足:

\[\begin{equation*} |x_n-x^{\star}|\leqslant\frac{b-a}{2^n} \end{equation*}\]

因此,可通过给定误差容许值\(\varepsilon\),选取满足:

\[\begin{equation*} \frac{b-a}{2^n}\leqslant\varepsilon \end{equation*}\]

的最小整数\(n\)来控制计算精度。

14.1.0.2 不动点迭代法

Method 14.2. (Fixed-point Iteration Method)

为求解非线性方程\(f(x)=0\),可将其改写为等价的形式\(x=T(x)\),若\(T\)是一个压缩映射,由定理 3.20可求得非线性方程的根。对于迭代次数问题,可用定理 3.20中的第二种误差估计公式来控制精度。

Definition 14.1. 设映射\(T\)有不动点\(x^{\star}\)。若存在\(x^{\star}\)的邻域\(\bar{U}(x^{\star},\delta)\)满足对任意的\(x\in\bar{U}(x^{\star},\delta)\),序列\(\{x_n=T^nx\}\to x^{\star}\),则称此时的迭代法局部收敛(local convergence)

Theorem 14.1. 设映射\(T\)有不动点\(x^{\star}\)\(T'\)\(x^{\star}\)的某个邻域内连续,且\(|T'x^{\star}|<1\),则此时的迭代法局部收敛。

证明. 由条件和性质 3.2.9(5)可知存在\(x^{\star}\)的邻域\(U\)\(L\),使得对任意的\(x\in U\)\(|T'x|\leqslant L<1\)。根据定理 4.3可知对任意的\(x'\in U\)有:

\[\begin{equation*} |Tx-Tx^{\star}|\leqslant L|x-x^{\star}| \end{equation*}\]

\(T\)\(x^{\star}\)的某个闭邻域内是压缩映射,由定理 3.8性质 3.2.6(1)定理 3.20可知此时的迭代法局部收敛。 ◻

Definition 14.2. 设迭代过程\(x_{n+1}=Tx_n\)收敛于\(T\)的不动点\(x^{\star}\)。令\(\varepsilon_n=x_n-x^{\star}\),若有:

\[\begin{equation*} \lim_{n\to+\infty}\frac{\varepsilon_{n+1}}{\varepsilon_n^p}=C\in\mathbb{R}^{}\setminus\{0\} \end{equation*}\]

则称该迭代过程是\(p\)阶收敛的。特别的,\(p=1\)\(|C|<1\)时称为线性收敛(linear convergence)\(p>1\)时称为超线性收敛(superlinear convergence)

Theorem 14.2. 设迭代过程\(x_{n+1}=Tx_n\)。若对于\(p\in\mathbb{N}^+\)\(T^{(p)}\)\(T\)的不动点\(x^{\star}\)的某个邻域内连续,并且有:

\[\begin{equation*} T'x^{\star}=T''x^{\star}=\cdots=T^{(p-1)}x^{\star}=0,\;T^{(p)}x^{\star}\ne0 \end{equation*}\]

则该迭代过程在\(x^{\star}\)邻近是\(p\)阶收敛的。

证明. \(T'x^{\star}=0\)定理 14.1可知该迭代过程是局部收敛的。根据泰勒展开Lagrange余项可得:

\[\begin{equation*} Tx_n=Tx^{\star}+\frac{T^{(p)}\xi_n}{p!}(x_n-x^{\star})^p,\quad\xi_n\text{在$x_n$与$x^{\star}$之间} \end{equation*}\]

所以:

\[\begin{equation*} \frac{Tx_n-x^{\star}}{(x_n-x^{\star})^p}=\frac{T^{(p)}\xi_n}{p!} \end{equation*}\]

\(T^{(p)}\)\(x^{\star}\)某个邻域内的连续性和性质 3.2.9(8.c)(4.a)可得:

\[\begin{equation*} \lim_{n\to+\infty}\frac{Tx_n-x^{\star}}{(x_n-x^{\star})^p}=\lim_{n\to+\infty}\frac{T^{(p)}\xi_n}{p!}=\frac{1}{p!}T^{(p)}\left(\lim_{n\to+\infty}\xi_n\right)=\frac{T^{(p)}x^{\star}}{p!} \end{equation*}\]

14.1.0.3 牛顿法

Method 14.3. (Newton Method)

\(f\)在根\(x^{\star}\)的某个邻域内是\(C^2\)类函数,且满足\(f(x^{\star})=0,\;f'(x^{\star})\ne0\)。为求解非线性方程\(f(x)=0\),牛顿法对\(f\)在当前迭代点\(x_n\)处作一阶Taylor展开,得到:

\[\begin{equation*} f(x)\approx f(x_n)+f'(x_n)(x-x_n) \end{equation*}\]

令上式为\(0\)得到迭代格式为:

\[\begin{equation*} x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)} \end{equation*}\]

牛顿法的迭代函数为:

\[\begin{equation*} Tx=x-\frac{f(x)}{f'(x)} \end{equation*}\]

于是有:

\[\begin{equation*} T'x=1-\frac{[f'(x)]^2-f(x)f''(x)}{[f'(x)]^2}=\frac{f(x)f''(x)}{[f'(x)]^2} \end{equation*}\]

根据定理 14.2可知牛顿法至少是二阶收敛的。

14.2 插值法

Definition 14.3. 设函数\(f\)\([a,b]\)上有定义,且已知在点\(a\leqslant x_0<x_1<\cdots<x_n\leqslant b\)上的值\(y_0,y_1,\dots,y_n\)。若存在一个简单函数\(P(x)\)使得:

\[\begin{equation*} P(x_i)=y_i,\quad i=0,1,\dots,n \end{equation*}\]

则称\(P\)\(f\)插值函数(interpolating function),点\(x_0,x_1,\dots,x_n\)被称为插值节点(interpolation nodes)\([a,b]\)插值区间(interpolation interval),求插值函数\(P\)的方法称为插值法(interpolation method)。若\(P\)是不超过\(n\)次的多项式,即:

\[\begin{equation*} P(x)=\sum_{i=0}^{n}a_ix^i,\quad a_i\in\mathbb{R}^{},\;i=0,1,\dots,n \end{equation*}\]

则称\(P\)插值多项式(interpolation polynomial),对应的插值法被称为多项式插值(polynomial interpolation)。若\(P\)是分段的多项式,则称对应的插值法为分段插值(piecewise interpolation)

Theorem 14.3. 设函数\(f\)\([a,b]\)上有定义,且已知在点\(a\leqslant x_0<x_1<\cdots<x_n\leqslant b\)上的值\(y_0,y_1,\dots,y_n\)\(f\)\([a,b]\)上次数不超过\(n\)的插值多项式存在且唯一。

证明. \(P(x)=\sum\limits_{i=0}^{n}a_ix^i\),则由条件可得到:

\[\begin{equation*} \begin{pmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n \\ 1 & x_1 & x_1^2 & \cdots & x_1^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{pmatrix} = \begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_n \end{pmatrix} \end{equation*}\]

根据性质 2.1.7(13)性质 1.1.3(1)定理 2.4(2)可知上述线性方程组的解存在且唯一,即\(f\)\([a,b]\)上次数不超过\(n\)的插值多项式存在且唯一。 ◻

14.2.0.1 Lagrange插值多项式

Definition 14.4. 若\(n\)次多项式\(l_i(x),\;i=0,1,\dots,n\)\(n+1\)个节点\(x_0<x_1<\cdots<x_n\)上满足:

\[\begin{equation*} l_i(x_j)= \begin{cases} 1,& i=j \\ 0,& i\ne j \end{cases},\quad i,j=0,1,\dots,n \end{equation*}\]

则称\(l_0(x),l_1(x),\dots,l_n(x)\)为节点\(x_0,x_1,\dots,x_n\)上的\(n\)插值基函数(interpolation basis functions)

Definition 14.5. 设函数\(f\)\([a,b]\)上有定义,且已知在点\(a\leqslant x_0<x_1<\cdots<x_n\leqslant b\)上的值\(y_0,y_1,\dots,y_n\)。定义:

\[\begin{equation*} l_i(x)=\frac{(x-x_0)(x-x_1)\dots(x-x_{i-1})(x-x_{i+1})\dots(x-x_n)}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)},\quad i=0,1,\dots,n \end{equation*}\]

称:

\[\begin{equation*} L_n(x)\coloneq\sum_{i=0}^{n}y_il_i(x) \end{equation*}\]

\(f\)在节点\(x_0,x_1,\dots,x_n\)上的拉格朗日插值多项式(Lagrange interpolation polynomial)

Property 14.2.1. 设函数\(f\)\([a,b]\)上有定义,且已知在点\(a\leqslant x_0<x_1<\cdots<x_n\leqslant b\)上的值\(y_0,y_1,\dots,y_n\)\(f\in C^n([a,b])\)\(f^{(n+1)}\)\((a,b)\)上存在,\(L_n\)\(f\)在节点\(x_0,x_1,\dots,x_n\)上的Lagrange插值多项式。\(L_n\)具有如下性质:

  1. \(w_{n+1}(x)=(x-x_0)(x-x_1)\cdots(x-x_n)\),则:

    \[\begin{equation*} L_n(x)=\sum_{i=0}^{n}y_i\frac{w_{n+1}(x)}{(x-x_i)w_{n+1}'(x_i)} \end{equation*}\]

  2. 对于任意的\(x\in[a,b]\),有误差估计:

    \[\begin{equation*} f(x)-L_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}w_{n+1}(x) \end{equation*}\]

    其中\(\xi\in(a,b)\)

  3. \(\sum\limits_{i=0}^{n}x_i^ml_i(x)=x^m(m\leqslant n)\)

  4. \(L_n\)有如下形式:

    \[\begin{equation*} L_n(x)=\frac{\sum\limits_{i=0}^{n}y_i\dfrac{1}{(x-x_i)w_{n+1}'(x_i)}}{\sum\limits_{i=0}^{n}\dfrac{1}{(x-x_i)w_{n+1}'(x_i)}} \end{equation*}\]

证明. (1)由乘积的微分可知:

\[\begin{equation*} w_{n+1}'(x_i)=(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n) \end{equation*}\]

于是:

\[\begin{equation*} l_i(x)=\frac{w_{n+1}(x)}{(x-x_i)w_{n+1}'(x_i)} \end{equation*}\]

(2)由条件可设\(f(x)-L_n(x)=K(x)w_{n+1}(x)\),其中\(K(x)\)是一个待定函数。作辅助函数:

\[\begin{equation*} \varphi(t)=f(t)-L_n(t)-K(x)w_{n+1}(t) \end{equation*}\]

性质 3.3.3(5.b)性质 4.1.1(4)可知\(\varphi^{(n)}\)\([a,b]\)上连续,\(\varphi^{(n+1)}\)\((a,b)\)上存在。因为\(\varphi\)\(x_0,x_1,\dots,x_n\)\(x\)处均为\(0\),所以\(\varphi\)\([a,b]\)上有\(n+2\)个零点,根据性质 3.3.3(5.b)性质 4.1.1(4)定理 4.2可知\(\varphi'\)\(\varphi\)的两个零点间至少存在一个零点,所以\(\varphi'\)\([a,b]\)上至少存在\(n+1\)个零点。对\(\varphi'\)再使用定理 4.2可知\(\varphi''\)\([a,b]\)上至少存在\(n\)个零点。依此类推可知\(\varphi^{(n+1)}\)\((a,b)\)上至少存在一个零点\(\xi\),即:

\[\begin{equation*} \varphi^{(n+1)}(\xi)=f^{(n+1)}(\xi)-K(x)(n+1)!=0 \end{equation*}\]

于是:

\[\begin{equation*} K(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!} \end{equation*}\]

即:

\[\begin{equation*} f(x)-L_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}w_{n+1}(x) \end{equation*}\]

(3)设\(f=x^m\),由(2)即可得出结论。

(4)在(3)中取\(m=0\),由(1)(3)可得:

\[\begin{equation*} 1=\sum_{i=0}^{n}l_i(x)=\sum_{i=0}^{n}\frac{w_{n+1}(x)}{(x-x_i)w_{n+1}'(x_i)}=w_{n+1}(x)\sum_{i=0}^{n}\frac{1}{(x-x_i)w_{n+1}'(x_i)} \end{equation*}\]

所以:

\[\begin{equation*} w_{n+1}(x)=\frac{1}{\sum\limits_{i=0}^{n}\dfrac{1}{(x-x_i)w_{n+1}'(x_i)}} \end{equation*}\]

于是有:

\[\begin{align*} L_n(x)&=\sum_{i=0}^{n}y_i\frac{w_{n+1}(x)}{(x-x_i)w_{n+1}'(x_i)}=w_{n+1}(x)\sum_{i=0}^{n}y_i\frac{1}{(x-x_i)w_{n+1}'(x_i)} \\ &=\frac{\sum\limits_{i=0}^{n}y_i\dfrac{1}{(x-x_i)w_{n+1}'(x_i)}}{\sum\limits_{i=0}^{n}\dfrac{1}{(x-x_i)w_{n+1}'(x_i)}} \end{align*}\]

note 14.1. 性质 14.2.1(4)\(L_n\)的形式又被称为Barycentric形式。该表示形式在数学上与拉格朗日插值多项式完全等价,但在数值计算中避免了高次多项式乘积的显式计算,即计算:

\[\begin{equation*} w_{n+1}(x)=(x-x_0)(x-x_1)\cdots(x-x_n) \end{equation*}\]

从而减少了由大数乘积和相消所引起的数值误差,具有显著更优的稳定性。基于上述原因,在实际数值计算中通常优先采用Barycentric形式来实现多项式插值。

    import numpy as np
    from scipy.interpolate import BarycentricInterpolator
    x_nodes = np.array([0.0,1.0,2.0,3.0])
    y_nodes = np.array([1.0,2.0,0.0,5.0])
    interp = BarycentricInterpolator(x_nodes,y_nodes)
    print(interp(1.5))

14.2.0.2 Newton插值多项式

Definition 14.6. 设函数\(f\)\([a,b]\)上有定义,\(a\leqslant x_0<x_1<\cdots<x_n\leqslant b\)。记\(f[x_0]=f(x_0)\),称:

\[\begin{equation*} f[x_0,x_n]=\frac{f(x_n)-f(x_0)}{x_n-x_0} \end{equation*}\]

\(f\)关于点\(x_0,x_n\)一阶均差(first-order divided difference)。称:

\[\begin{equation*} f[x_0,x_1,\dots,x_n]=\frac{f[x_1,x_2,\dots,x_n]-f[x_0,x_1,\dots,x_{n-1}]}{x_n-x_0} \end{equation*}\]

\(f\)n阶均差(n-th order divided difference)

Property 14.2.2. 设函数\(f\)\([a,b]\)上有定义,\(a\leqslant x_0<x_1<\cdots<x_n\leqslant b\)。均差具有如下性质:

  1. \(f\)\(n\)阶均差可以表示为\(f(x_0),f(x_1),\dots,f(x_n)\)的线性组合:

    \[\begin{equation*} f[x_0,x_1,\dots,x_n]=\sum_{i=0}^{n}\frac{f(x_i)}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)} \end{equation*}\]

  2. \(f[x_0,x_1,\dots,x_n]\)\(\{x_0,x_1,\dots,x_n\}\)的任意顺序置换保持不变;

  3. 对任意的\(\alpha,\beta\in\mathbb{R}^{}\)\([a,b]\)上的函数\(g\),有:

    \[\begin{equation*} (\alpha f+\beta g)[x_0,x_1,\dots,x_n]=\alpha f[x_0,x_1,\dots,x_n]+\beta g[x_0,x_1,\dots,x_n] \end{equation*}\]

  4. \(f\in C^n([a,b])\),则存在\(\xi\in[a,b]\)使得:

    \[\begin{equation*} f[x_0,x_1,\dots,x_n]=\frac{f^{(n)}(\xi)}{n!} \end{equation*}\]

证明. (1)使用数学归纳法。当\(n=1\)时有:

\[\begin{equation*} f[x_0,x_1]=\frac{f(x_1)-f(x_0)}{x_1-x_0}=\frac{f(x_0)}{x_0-x_1}+\frac{f(x_1)}{x_1-x_0} \end{equation*}\]

结论成立。设结论对\(n-1\)成立,下面证明结论对\(n\)也成立。

根据归纳假设:

\[\begin{align*} &f[x_0,x_1,\dots,x_n]=\frac{f[x_1,x_2,\dots,x_n]-f[x_0,x_1,\dots,x_{n-1}]}{x_n-x_0} \\ =&\frac{1}{x_n-x_0}\left[\sum_{i=1}^{n}\frac{f(x_i)}{(x_i-x_1)(x_i-x_2)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)}\right. \\ &\left.-\sum_{i=0}^{n-1}\frac{f(x_i)}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_{n-1})}\right] \\ =&\frac{f(x_n)}{(x_n-x_0)(x_n-x_1)\dots(x_n-x_{n-1})}+\frac{f(x_0)}{(x_0-x_1)(x_0-x_2)\dots(x_0-x_n)} \\ &+\sum_{i=1}^{n-1}\left[\frac{f(x_i)}{(x_n-x_0)(x_i-x_1)(x_i-x_2)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)}\right. \\ &\left.-\frac{f(x_i)}{(x_n-x_0)(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_{n-1})}\right] \end{align*}\]

对于最后一项:

\[\begin{align*} &\sum_{i=1}^{n-1}\left[\frac{(x_i-x_0)f(x_i)}{(x_n-x_0)(x_i-x_0)(x_i-x_1)(x_i-x_2)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)}\right. \\ &\left.-\frac{(x_i-x_n)f(x_i)}{(x_n-x_0)(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_{n-1})(x_i-x_n)}\right] \\ =&\sum_{i=1}^{n-1}\frac{(x_i-x_0)f(x_i)-(x_i-x_n)f(x_i)}{(x_n-x_0)(x_i-x_0)(x_i-x_1)(x_i-x_2)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)} \\ =&\sum_{i=1}^{n-1}\frac{(x_n-x_0)f(x_i)}{(x_n-x_0)(x_i-x_0)(x_i-x_1)(x_i-x_2)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)} \\ =&\sum_{i=1}^{n-1}\frac{f(x_i)}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)} \end{align*}\]

所以结论成立。

(2)由(1)立即可得。

(3)由(1)立即可得。

(4)令\(L_{n-1}\)\(f\)在节点\(x_0,x_1,\dots,x_{n-1}\)上的Lagrange插值多项式,作辅助函数:

\[\begin{equation*} \varphi(x)=f(x)-L_{n-1}(x)-f[x_0,x_1,\dots,x_n]w_n(x) \end{equation*}\]

\(L_{n-1}\)的定义可知\(f-L_{n-1}\)在节点\(x_0,x_1,\dots,x_{n-1}\)上为\(0\),而\(w_n\)在节点\(x_0,x_1,\dots,x_{n-1}\)上为\(0\),所以\(x_0,x_1,\dots,x_{n-1}\)\(\varphi\)的零点。注意到:

\[\begin{align*} &f(x_n)-L_{n-1}(x_n)-f[x_0,x_1,\dots,x_n]w_n(x_n) \\ =&f(x_n)-\sum_{i=0}^{n-1}\frac{f(x_i)(x_n-x_0)(x_n-x_1)\dots(x_n-x_{i-1})(x_n-x_{i+1})\dots(x_n-x_{n-1})}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_{n-1})} \\ &-\sum_{i=0}^{n}\frac{f(x_i)(x_n-x_0)(x_n-x_1)\cdots(x_n-x_{n-1})}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)} \\ =&f(x_n)-\sum_{i=0}^{n-1}\frac{f(x_i)(x_n-x_0)(x_n-x_1)\dots(x_n-x_{i-1})(x_n-x_{i+1})\dots(x_n-x_{n-1})}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_{n-1})} \\ &-f(x_n)-\sum_{i=0}^{n-1}\frac{f(x_i)(x_n-x_0)(x_n-x_1)\cdots(x_n-x_{n-1})}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)} \\ =&-\sum_{i=0}^{n-1}\frac{f(x_i)(x_n-x_0)(x_n-x_1)\dots(x_n-x_{n-1})}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_{n-1})(x_n-x_i)} \\ &-\sum_{i=0}^{n-1}\frac{f(x_i)(x_n-x_0)(x_n-x_1)\cdots(x_n-x_{n-1})}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)}=0 \end{align*}\]

所以\(x_n\)也是\(\varphi\)的零点。根据性质 3.3.3(5.b)性质 4.1.1(4)定理 4.2可知\(\varphi'\)\([a,b]\)上存在\(n\)个零点。对\(\varphi'\)再使用定理 4.2可知\(\varphi''\)\([a,b]\)上存在\(n-1\)个零点。依此类推可知\(\varphi^{(n)}\)\([a,b]\)上存在一个零点\(\xi\),即:

\[\begin{equation*} \varphi^{(n)}(\xi)=f^{(n)}(\xi)-f[x_0,x_1,\dots,x_n]n!=0 \end{equation*}\]

所以:

\[\begin{equation*} f[x_0,x_1,\dots,x_n]=\frac{f^{(n)}(\xi)}{n!} \end{equation*}\]

Definition 14.7. 设函数\(f\)\([a,b]\)上有定义,且已知在点\(a\leqslant x_0<x_1<\cdots<x_n\leqslant b\)上的值\(y_0,y_1,\dots,y_n\)。称:

\[\begin{equation*} P_n(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)+\cdots+f[x_0,x_1,\dots,x_n]\prod_{i=0}^{n-1}(x-x_i) \end{equation*}\]

\(f\)在节点\(x_0,x_1,\dots,x_n\)上的牛顿插值多项式(Newton interpolation polynomial)

Property 14.2.3. Newton插值多项式与Lagrange插值多项式等价。

证明. 设函数\(f\)\([a,b]\)上有定义,且已知在点\(a\leqslant x_0<x_1<\cdots<x_n\leqslant b\)上的值\(y_0,y_1,\dots,y_n\)。因为\(f\)在节点\(x_0,x_1,\dots,x_n\)上的Newton插值多项式\(P_n\)的次数不超过\(n\),根据定理 14.3只需证明\(P_n(x_i)=f(x_i),\;i=0,1,\dots,n\)

使用数学归纳法。当\(i=0\)时由\(P_n\)的定义可知\(P_n(x_0)=f(x_0)\),设结论对\(n-1\)成立,下面证明结论对\(n\)也成立。

对于:

\[\begin{equation*} P_n(x)=P_{n-1}(x)+f[x_0,x_1,\dots,x_n]w_n(x) \end{equation*}\]

因为对任意的\(i=0,1,\dots,n-1\)\(w_n(x_i)=0\),所以由归纳假设可知:

\[\begin{equation*} P_n(x_i)=P_{n-1}(x_i)=f(x_i),\quad i=0,1,\dots,n-1 \end{equation*}\]

下面证明\(P_n(x_n)=f(x_n)\)

作辅助函数\(R(x)=f(x)-P_{n-1}(x)\),于是有:

\[\begin{equation*} \frac{R(x_n)}{w_n(x_n)}=\frac{f(x_n)-P_{n-1}(x_n)}{w_n(x_n)} \end{equation*}\]

性质 14.2.2(1)可得:

\[\begin{equation*} R[x_0,x_1,\dots,x_n]=\sum_{i=0}^{n}\frac{R(x_i)}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)}=\frac{R(x_n)}{w_n(x_n)} \end{equation*}\]

所以:

\[\begin{equation*} \frac{f(x_n)-P_{n-1}(x_n)}{w_n(x_n)}=R[x_0,x_1,\dots,x_n] \end{equation*}\]

根据性质 14.2.2(3)可知:

\[\begin{equation*} R[x_0,x_1,\dots,x_n]=f[x_0,x_1,\dots,x_n]-P_{n-1}[x_0,x_1,\dots,x_n] \end{equation*}\]

因为\(P_{n-1}\)的次数不超过\(n-1\),由性质 14.2.2(4)可知\(P_{n-1}[x_0,x_1,\dots,x_n]=0\),于是有:

\[\begin{equation*} R[x_0,x_1,\dots,x_n]=f[x_0,x_1,\dots,x_n] \end{equation*}\]

所以:

\[\begin{equation*} \frac{f(x_n)-P_{n-1}(x_n)}{w_n(x_n)}=f[x_0,x_1,\dots,x_n] \end{equation*}\]

即:

\[\begin{equation*} P_{n-1}(x_n)=f(x_n)-f[x_0,x_1,\dots,x_n]w_n(x_n) \end{equation*}\]

于是:

\[\begin{equation*} P_n(x_n)=P_{n-1}(x_n)+f[x_0,x_1,\dots,x_n]w_n(x_n)=f(x_n) \end{equation*}\]

note 14.2. 根据性质 14.2.3,Newton插值多项式和Lagrange插值多项式是等价的,于是它具有了Lagrange插值多项式的全部性质。既然等价,那为什么要研究它呢?

从构造方式看,Lagrange插值是一次性使用全部节点进行整体构造,而Newton插值具有天然的递推性。当新增插值节点时,Lagrange插值多项式通常需要重新计算,而Newton插值多项式只需增加一个新的均差项即可完成更新,这使其在逐点插值与在线计算中具有明显优势。

14.2.0.3 Hermite插值

Definition 14.8. 设函数\(f\)\([a,b]\)上有定义,且已知\(f^{(j)}(x_i)\)\(i=0,1,\dots,n,\;j=0,1,\dots,m_i\),其中\(a\leqslant x_0<x_1<\cdots<x_n\leqslant b\)。令\(N=\sum\limits_{i=0}^{n}m_i-1\),称满足:

\[\begin{equation*} H^{(j)}_N(x_i)=f^{(j)}(x_i),\quad i=0,1,\dots,n,\;j=0,1,\dots,m_i \end{equation*}\]

的不超过\(N\)次的多项式\(H_N(x)\)埃尔米特插值多项式(Hermite interpolation polynomial)

Definition 14.9. 定义:

\[\begin{equation*} f[\underbrace{x,x,\dots,x}_{n+1\text{个}x}]\coloneq\frac{f^{(n)}(x)}{n!} \end{equation*}\]

Property 14.2.4. 设函数\(f\)\([a,b]\)上有定义,且已知\(f^{(j)}(x_i)\)\(i=0,1,\dots,n,\;j=0,1,\dots,m_i\),其中\(a\leqslant x_0<x_1<\cdots<x_n\leqslant b\)\(N=\sum\limits_{i=0}^{n}m_i-1\)。Hermite插值多项式具有如下性质:

  1. Hermite插值多项式存在且唯一;

  2. 若将每个节点\(x_i\)重复\(m_i+1\)次,按重复顺序排列为\(y_0,y_1,\dots,y_M\),则:

    \[\begin{equation*} H_N(x)=f[y_0]+f[y_0,y_1](x-y_0)+\cdots+f[y_0,y_1,\dots,y_M]\prod_{i=0}^{M-1}(x-y_i) \end{equation*}\]

    import numpy as np
    from scipy.interpolate import KroghInterpolator
    # Hermite interpolation data
    # f(0) = 1, f'(0) = 0
    # f(1) = 2, f'(1) = -1
    x = np.array([0.0, 0.0, 1.0, 1.0])
    y = np.array([1.0, 0.0, 2.0, -1.0])
    # Hermite interpolation
    hermite_poly = KroghInterpolator(x, y)
    # evaluate
    x_eval = np.linspace(0.0, 1.0, 100)
    y_eval = hermite_poly(x_eval)
    print(hermite_poly(0.5))

14.3 变分法

Definition 14.10. 设\(\varphi,\mu\)是可测空间\((X,\mathscr{F})\)上的概率测度,\(\varphi\ll\mu\),根据定理 5.35,称:

\[\begin{equation*} \int_{X}\ln\frac{\mathop{}\!\mathrm{d}\varphi}{\mathop{}\!\mathrm{d}\mu}(x)\mathop{}\!\mathrm{d}\varphi \end{equation*}\]

\(\varphi\)关于\(\mu\)KL散度(Kullback–Leibler divergence) ,记为\(\operatorname{KL}(\varphi\|\mu)\)

note 14.3. 设\(\mu\)是可测空间\((X,\mathscr{F})\)上的概率测度,我们希望求解\(\mu\)但是它很难求解,变分法(variational method)的思想就是取一族形式简单、便于计算的概率函数,并在这一族分布中找到一个最接近\(p\)的。由不等式 20可以看出KL散度可以被用来衡量分布之间的差异性。

note 14.4. 在贝叶斯统计中,我们需要求解后验分布\(p(y\mid x)\),变分法提供了一个视角。

\(y\mid x\)的概率分布为\(\varphi\),其近似概率分布为\(\psi\)\(\varphi,\psi\ll\mu\)\(\varphi,\psi,\mu\)是可测空间\((X,\mathscr{F})\)上的概率测度,根据定理 5.35,记:

\[\begin{equation*} q(y)=\frac{\mathop{}\!\mathrm{d}\psi}{\mathop{}\!\mathrm{d}\mu},\quad p(y\mid x)=\frac{\mathop{}\!\mathrm{d}\varphi}{\mathop{}\!\mathrm{d}\mu} \end{equation*}\]

性质 5.6.3(5)引理 5.10可知:

\[\begin{equation*} \operatorname{KL}(\psi\|\varphi)=\int_{X}\ln\frac{\mathop{}\!\mathrm{d}\psi}{\mathop{}\!\mathrm{d}\varphi}\mathop{}\!\mathrm{d}\psi=\int_{X}\ln\left(\frac{\mathop{}\!\mathrm{d}\psi}{\mathop{}\!\mathrm{d}\mu}\Big/\frac{\mathop{}\!\mathrm{d}\varphi}{\mathop{}\!\mathrm{d}\mu}\right)\mathop{}\!\mathrm{d}\psi=\int_{X}\ln\frac{q(y)}{p(y\mid x)}q(y)\mathop{}\!\mathrm{d}\mu \end{equation*}\]

由于KL散度公式中存在未知的\(p(y\mid x)\),于是需要将其转变为易于求解的形式。由贝叶斯公式性质 5.4.3(6)可得:

\[\begin{align*} &\operatorname{KL}(\psi\|\varphi)=\int_{X}q(y)\ln\frac{q(y)}{p(y\mid x)}\mathop{}\!\mathrm{d}\mu=\int_{X}q(y)[\ln q(y)-\ln p(y\mid x)]\mathop{}\!\mathrm{d}\mu \\ =&\int_{X}q(y)[\ln q(y)-\ln p(y,x)+\ln p(x)]\mathop{}\!\mathrm{d}\mu \\ =&\int_{X}q(y)[\ln q(y)-\ln p(y,x)]\mathop{}\!\mathrm{d}\mu+\int_{X}q(y)\ln p(x)\mathop{}\!\mathrm{d}\mu \\ =&\int_{X}q(y)[\ln q(y)-\ln p(y,x)]\mathop{}\!\mathrm{d}\mu+\ln p(x)\int_{X}q(y)\mathop{}\!\mathrm{d}\mu \\ =&\int_{X}q(y)[\ln q(y)-\ln p(y,x)]\mathop{}\!\mathrm{d}\mu+\ln p(x) \end{align*}\]

于是根据性质 5.4.3(6)可得:

\[\begin{equation*} \ln p(x)=\int_{X}q(y)[\ln p(y,x)-\ln q(y)]\mathop{}\!\mathrm{d}\mu+\operatorname{KL}(\psi\|\varphi) \end{equation*}\]

因为\(\ln p(x)\)\(y\)无关,所以最小化\(\operatorname{KL}(q\| p)\)等价于最大化\(\int_{X}q(y)[\ln p(y,x)-\ln q(y)]\mathop{}\!\mathrm{d}\mu\)

Definition 14.11. 在前述框架下,称:

\[\begin{equation*} \int_{X}q(y)[\ln p(y,x)-\ln q(y)]\mathop{}\!\mathrm{d}\mu \end{equation*}\]

证据下界(evidence lower bound),记做\(L(q,x)\)

14.3.1 EM算法

在极大似然估计中,若观测变量\(f\)为从可控参数结构\((X,\mathscr{A},\mathscr{P})\)到可测空间\((Y,\mathscr{B})\)上的可测映射,参数空间为\(\Theta\),控制测度为\(\mu\),根据定理 5.35,记\(\dfrac{\mathop{}\!\mathrm{d}P_{\theta}}{\mathop{}\!\mathrm{d}\mu}(f)=p_{\theta}(f)\),那么目标函数可以写成:

\[\begin{equation*} \hat{\theta}=\underset{\theta}{\arg\max}\ln p_{\theta}(f) \end{equation*}\]

然而,在许多含有潜在结构的统计模型中,观测变量\(f\)本身的概率分布并不容易直接写出。真正容易描述的,往往是它在某些无法观测的隐变量\(g\)给定条件下的条件分布。有时,隐变量\(g\)本身的分布也不是单纯的辅助工具,而是我们希望估计和理解的对象。

但是在引入隐变量之后,我们需要通过对隐变量积分,才能得到观测变量的边际似然。设\(g\)为概率空间\((X,\mathscr{A},P)\)到可测空间\((Z,\mathscr{C})\)上的可测映射,\(Pg^{-1}\ll\nu\),根据定理 5.35,记\(\dfrac{\mathop{}\!\mathrm{d}Pg^{-1}}{\mathop{}\!\mathrm{d}\nu}=q(g)\),由全概率公式引理 5.10可得:

\[\begin{equation*} p_{\theta}(f)=\int_{Z}p_{\theta}(f\mid g)\mathop{}\!\mathrm{d}Pg^{-1}=\int_{Z}q(g)p_{\theta}(f\mid g)\mathop{}\!\mathrm{d}\nu \end{equation*}\]

原本我们只是想最大化一个对数似然函数,但由于隐变量无法观测,观测数据的概率必须通过积分得到,这使得目标函数通常难以直接优化。EM算法(Expectation-Maximization Algorithm)正是为了解决这类问题而提出的。

不等式 19贝叶斯公式性质 5.4.3(6)可得:

\[\begin{align*} &\ln p_{\theta}(f)=\ln\int_{Z}q(g)p_{\theta}(f\mid g)\mathop{}\!\mathrm{d}\nu=\ln\operatorname{E}[p_{\theta}(f\mid g)]\geqslant\operatorname{E}[\ln p_{\theta}(f\mid g)] \\ =&\int_{Z}q(g)[\ln p_{\theta}(f\mid g)]\mathop{}\!\mathrm{d}\nu=\int_{Z}q(g)[\ln p_{\theta}(f,g)-\ln q(g)]\mathop{}\!\mathrm{d}\nu \\ =&\int_{Z}q(g)\ln p_{\theta}(f,g)\mathop{}\!\mathrm{d}\nu-\int_{Z}q(g)\ln q(g)\mathop{}\!\mathrm{d}\nu \end{align*}\]

记:

\[\begin{equation*} L(q,y)=\int_{Z}q(g)\ln p_{\theta}(f,g)\mathop{}\!\mathrm{d}\nu-\int_{Z}q(g)\ln q(g)\mathop{}\!\mathrm{d}\nu \end{equation*}\]

则有:

\[\begin{equation*} L(q,y)\leqslant\ln p_{\theta}(f) \end{equation*}\]

这个\(L(q,y)\)就是边际对数似然\(\ln p_{\theta}(f)\)的一个下界。

EM算法的基本思想就是:既然直接最大化\(\ln p_{\theta}(f)\)比较困难,那么可以转而反复提高这个下界。若一组\(\hat{q},\hat{\theta}\)使得\(L(q,y)\)达到最大,则\(\hat{\theta}\)也应该是\(\theta\)的极大似然估计。

进一步地,由性质 5.6.3(5)引理 5.10贝叶斯公式可知这个下界与真实边际对数似然之间的差距可以精确写成一个KL散度:

\[\begin{align*} &\operatorname{KL}[q(g)\|p_{\theta}(g\mid f)]=\int_{Z}q(g)\ln\frac{q(g)}{p_{\theta}(g\mid f)}\mathop{}\!\mathrm{d}\nu \\ =&\int_{Z}q(g)\ln q(g)\mathop{}\!\mathrm{d}\nu- \int_{Z}q(g)\ln p_{\theta}(g\mid f)\mathop{}\!\mathrm{d}\nu \\ =&\int_{Z}q(g)\ln q(g)\mathop{}\!\mathrm{d}\nu-\int_{Z}q(g)\ln p_{\theta}(f,g)\mathop{}\!\mathrm{d}\nu+\ln p_{\theta}(f) \end{align*}\]

于是:

\[\begin{equation*} \ln p_{\theta}(f)=L(q,y)+\operatorname{KL}[q(g)\|p_{\theta}(g\mid f)] \end{equation*}\]

Method 14.4. (Expectation-Maximization Algorithm)

EM算法可以理解为对\(L(q,y)\)进行坐标上升。具体地说,在第\(t\)轮迭代中,假设当前参数为\(\theta^{(t)}\)。EM算法交替更新\(q\)\(\theta\)

首先固定\(\theta=\theta^{(t)}\),求解\(q^{(t+1)}\)使得最大化\(L(q^{(t+1)},y)\)。由于:

\[\begin{equation*} \ln p_{\theta^{(t)}}(f)=L(q^{(t+1)},y)+\operatorname{KL}[q^{(t+1)}(g)\|p_{\theta^{(t)}}(g\mid f)] \end{equation*}\]

而左边与\(q\)无关,是一个定值,所以最大化\(L(q^{(t+1)},y)\)等价于最小化\(\operatorname{KL}[q^{(t+1)}(g)\|p_{\theta^{(t)}}(g\mid f)]\)。当\(q^{(t+1)}(g)=p_{\theta^{(t)}}(g\mid f)\)时,KL散度达到最小值\(0\)。因此,EM算法的E步就是在当前参数\(\theta^{(t)}\)下计算隐变量的后验分布。

接着固定\(q=q^{(t+1)}\),最大化\(L(q^{(t+1)},y)\)。因为:

\[\begin{equation*} L[q^{(t+1)},y]=\int_{Z}q^{(t+1)}(g)\ln p_{\theta}(f,g)\mathop{}\!\mathrm{d}\nu-\int_{Z}q^{(t+1)}(g)\ln q^{(t+1)}(g)\mathop{}\!\mathrm{d}\nu \end{equation*}\]

其中第二项不含\(\theta\),所以关于\(\theta\)的优化等价于:

\[\begin{equation*} \theta^{(t+1)}=\underset{\theta}{\arg\max}\int_{Z}q^{(t+1)}(g)\ln p_{\theta}(f,g)\mathop{}\!\mathrm{d}\nu \end{equation*}\]

代入:

\[\begin{equation*} q^{(t+1)}(g)=p_{\theta^{(t)}}(g\mid f) \end{equation*}\]

得到:

\[\begin{equation*} \theta^{(t+1)}=\underset{\theta}{\arg\max}\int_{Z}p_{\theta^{(t)}}(g\mid f)\ln p_{\theta}(f,g)\mathop{}\!\mathrm{d}\nu \end{equation*}\]

对于这个优化问题,我们就可以使用一般的极大似然估计去做了。

上述第一步被称之为\(E\)步,第二步被称之为\(M\)步,称:

\[\begin{equation*} \int_{Z}q^{(t+1)}(g)\ln p_{\theta}(f,g)\mathop{}\!\mathrm{d}\nu \end{equation*}\]

为Q函数。

COMPUTATIONAL NOTES

计算实验

实验01

非线性方程求解:迭代误差

Python · Jupyter

非线性方程:误差随迭代次数的变化

求解方程

cos x − x = 0.

使用 SciPy 的 Brent 法计算高精度参考根,再比较二分法、不动点迭代与 Newton 法的逐步误差。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import root_scalar

def f(x):
    return np.cos(x) - x

def df(x):
    return -np.sin(x) - 1.0

reference = root_scalar(
    f,
    bracket=(0.0, 1.0),
    method="brentq",
    xtol=1e-14,
).root
print(f"Reference root x* = {reference:.16f}")
Reference root x* = 0.7390851332151612
def bisection_history(left, right, iterations):
    values = []
    for _ in range(iterations):
        midpoint = (left + right) / 2.0
        values.append(midpoint)
        if f(left) * f(midpoint) <= 0:
            right = midpoint
        else:
            left = midpoint
    return np.asarray(values)

def fixed_point_history(x0, iterations):
    values = [x0]
    for _ in range(iterations):
        values.append(np.cos(values[-1]))
    return np.asarray(values)

def newton_history(x0, iterations):
    values = [x0]
    for _ in range(iterations):
        current = values[-1]
        values.append(current - f(current) / df(current))
    return np.asarray(values)

histories = {
    "Bisection": bisection_history(0.0, 1.0, 45),
    "Fixed-point iteration": fixed_point_history(0.5, 45),
    "Newton's method": newton_history(0.5, 8),
}

tiny = np.finfo(float).tiny
errors = {
    name: np.maximum(np.abs(values - reference), tiny)
    for name, values in histories.items()
}
colors = {
    "Bisection": "#52649a",
    "Fixed-point iteration": "#c26b4a",
    "Newton's method": "#397d78",
}

fig, ax = plt.subplots(figsize=(9.2, 5.4), constrained_layout=True)
for name, values in errors.items():
    ax.semilogy(
        np.arange(len(values)),
        values,
        marker="o",
        markersize=3.5,
        linewidth=1.8,
        color=colors[name],
        label=name,
    )

ax.set(
    title="Error in Nonlinear Equation Solvers",
    xlabel="Iteration $k$",
    ylabel=r"Absolute error $|x_k-x^*|$",
)
ax.grid(True, which="both", alpha=0.22)
ax.legend()
plt.show()

print("Method                    Iterates    Final absolute error")
for name, values in errors.items():
    print(f"{name:<12} {len(values):>6d}    {values[-1]:.3e}")

Method                    Iterates    Final absolute error
Bisection        45    2.276e-14
Fixed-point iteration     46    4.021e-09
Newton's method      9    5.551e-16

观察

  • 二分法误差稳定地按固定比例下降。
  • 不动点迭代为线性收敛,速度取决于根附近的压缩系数。
  • Newton 法进入根的邻域后误差近似平方下降,很快达到浮点精度。
实验02

多项式逼近:节点选择与逼近误差

Python · Jupyter

多项式插值:节点选择与逼近误差

在区间 [−1, 1] 上逼近 Runge 函数

$$f(x)=\frac{1}{1+25x^2}.$$

使用 SciPy 的 BarycentricInterpolator,比较等距节点与 Chebyshev 节点。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import BarycentricInterpolator

def runge(x):
    return 1.0 / (1.0 + 25.0 * x**2)

def equidistant_nodes(degree):
    return np.linspace(-1.0, 1.0, degree + 1)

def chebyshev_nodes(degree):
    k = np.arange(degree + 1)
    return np.cos((2 * k + 1) * np.pi / (2 * (degree + 1)))

grid = np.linspace(-1.0, 1.0, 3000)
truth = runge(grid)
degrees = [4, 8, 12, 16]
records = []

for degree in degrees:
    for node_name, node_function in [
        ("Equidistant nodes", equidistant_nodes),
        ("Chebyshev nodes", chebyshev_nodes),
    ]:
        nodes = node_function(degree)
        interpolator = BarycentricInterpolator(nodes, runge(nodes))
        approximation = interpolator(grid)
        maximum_error = np.max(np.abs(approximation - truth))
        records.append((degree, node_name, maximum_error))

print("Degree    Node type                 Maximum absolute error")
for degree, node_name, maximum_error in records:
    print(f"{degree:>4d}    {node_name:<18} {maximum_error:.4e}")
Degree    Node type                 Maximum absolute error
   4    Equidistant nodes  4.3836e-01
   4    Chebyshev nodes    4.0202e-01
   8    Equidistant nodes  1.0452e+00
   8    Chebyshev nodes    1.7084e-01
  12    Equidistant nodes  3.6634e+00
  12    Chebyshev nodes    6.9215e-02
  16    Equidistant nodes  1.4393e+01
  16    Chebyshev nodes    3.2613e-02
degree = 12
nodes_equal = equidistant_nodes(degree)
nodes_chebyshev = chebyshev_nodes(degree)

polynomial_equal = BarycentricInterpolator(
    nodes_equal,
    runge(nodes_equal),
)(grid)
polynomial_chebyshev = BarycentricInterpolator(
    nodes_chebyshev,
    runge(nodes_chebyshev),
)(grid)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), constrained_layout=True)

axes[0].plot(grid, truth, color="#1f2430", linewidth=2.2, label="Original function")
axes[0].plot(grid, polynomial_equal, color="#c26b4a", linewidth=1.7, label="Equidistant nodes")
axes[0].plot(grid, polynomial_chebyshev, color="#397d78", linewidth=1.7, label="Chebyshev nodes")
axes[0].scatter(nodes_equal, runge(nodes_equal), color="#c26b4a", s=18, zorder=3)
axes[0].set(
    title=f"Degree-{degree} Interpolation",
    xlabel="x",
    ylabel="Function value",
    ylim=(-0.5, 1.5),
)
axes[0].grid(alpha=0.2)
axes[0].legend()

axes[1].semilogy(
    grid,
    np.maximum(np.abs(polynomial_equal - truth), np.finfo(float).tiny),
    color="#c26b4a",
    linewidth=1.7,
    label="Equidistant nodes",
)
axes[1].semilogy(
    grid,
    np.maximum(np.abs(polynomial_chebyshev - truth), np.finfo(float).tiny),
    color="#397d78",
    linewidth=1.7,
    label="Chebyshev nodes",
)
axes[1].set(
    title="Pointwise Absolute Error",
    xlabel="x",
    ylabel="Absolute error",
)
axes[1].grid(True, which="both", alpha=0.2)
axes[1].legend()
plt.show()

观察

提高多项式次数并不必然改善全区间逼近。等距节点在区间端点附近可能出现明显振荡;Chebyshev 节点在端点附近更密集,通常能显著降低最大误差。