数值分析-函数逼近计算

前言

​   上一篇在多项式插值中我们介绍了几种经典方法，如Lagrange、Neville、Newton、逐次线性和Hermite插值等。多项式插值虽然可以在插值节点处一定经过给定的函数值，但是在插值节点之外的误差无法保证，经常出现剧烈震荡的情况。本章函数逼近计算的主要目的是在某种误差度量下求得最佳的逼近函数。多项式插值可以当作是一种局部准确拟合，函数逼近更多的是要求整体的拟合误差要小。

​   具体内容详见 《数值分析》 李庆扬、王能超、易大义（华中科大出版社）第三章

敲黑板

• 函数逼近其实很常见，在很多科学实验中的最小二乘拟合函数便是一个应用，在之前的矩阵论专题中有稍作讲解，相信大家都有一个目标：数据点都给我划到线上，啊啊啊！！

• 两种误差度量形式：设原函数为$f(x)$，x的取值范围为$a<x<b$，逼近函数为$P(x)$

  1. 一致逼近/均匀逼近：$||f(x)-P(x)|{|}_{\infty }=ma{x}_{a<x<b}|f(x)-P(x)|$
  2. 均方逼近/平方逼近：$||f(x)-P(x)|{|}_2=\sqrt{{\int }_a^b[f(x)-P(x){]}^{2}dx}$

• 根据教材的内容，我们来引入函数的内积空间及范数定义：

  • 函数的内积：设$f(x),g(x),\rho (x)\in C[a,b]$，其中$\rho (x)$是权函数

    权函数的定义为：${\int }_a^b|x{|}^n\rho (x)dx(n=0,1,..,n)$存在，且当${\int }_a^{b}g(x)\rho (x)dx=0$时，则$g(x)\equiv 0\in [a,b]$，其各阶矩也存在

    两个函数在区间$[a,b]$上的的内积表示为：$(f,g)={\int }_a^b\rho (x)f(x)g(x)dx$

    关于向量内积空间的性质同样适用于函数内积空间，如正定性、线性性、交换性

  • 函数的Euclid范数：$||f|{|}_2=\sqrt{{\int }_a^b\rho (x)f^2(x)dx}=\sqrt{(f,f)}$

    关于向量范数的不等式同样适用于函数范数，如柯西不等式、三角不等式、平行四边形定律。

  • 函数线性无关性：设${\varphi }_0(x),{\varphi }_1(x),...,{\varphi }_{n-1}(x)\in C[a,b]$，如果有

    $a_0{\varphi }_{(}x)+a_1{\varphi }_1(x)+...+a_{n-1}{\varphi }_{n-1}(x)=0$，当且仅当$a_0=a_1=...=a_{n-1}=0$时成立，则称${\varphi }_0(x),{\varphi }_1(x),...,{\varphi }_{n-1}(x)$在区间$[a,b]$上是线性无关的。

    给出判断函数线性相关性质的定理：${\varphi }_0(x),{\varphi }_1(x),...,{\varphi }_{n-1}(x)$在区间$[a,b]$上线性无关的充要条件是其克莱姆矩阵行列式不为0，即

    $$G_{n-1}=G({\varphi }_0(x),{\varphi }_1(x),...,{\varphi }_{n-1}(x))=\left|\begin{array}{cccc}({\varphi }_0,{\varphi }_0)& ({\varphi }_0,{\varphi }_1)& ...& ({\varphi }_0,{\varphi }_{n-1})\\ ({\varphi }_1,{\varphi }_0)& ({\varphi }_1,{\varphi }_1)& ...& ({\varphi }_1,{\varphi }_{n-1})\\ ⋮& ⋮& \ddots & ⋮\\ ({\varphi }_{n-1},{\varphi }_0)& ({\varphi }_{n-1},{\varphi }_1)& ...& ({\varphi }_{n-1},{\varphi }_{n-1})\end{array}\right|\text{\hspace{1em}(1)}$$

    • 这里试着证明其充分性，采用反证法，即若$G_{n-1}=0$，可以推得${\varphi }_0(x),{\varphi }_1(x),...,{\varphi }_{n-1}(x)$在区间$[a,b]$上线性相关：
      令
      $$A=\left[\begin{array}{cccc}({\varphi }_0,{\varphi }_0)& ({\varphi }_0,{\varphi }_1)& ...& ({\varphi }_0,{\varphi }_{n-1})\\ ({\varphi }_1,{\varphi }_0)& ({\varphi }_1,{\varphi }_1)& ...& ({\varphi }_1,{\varphi }_{n-1})\\ ⋮& ⋮& \ddots & ⋮\\ ({\varphi }_{n-1},{\varphi }_0)& ({\varphi }_{n-1},{\varphi }_1)& ...& ({\varphi }_{n-1},{\varphi }_{n-1})\end{array}\right]\text{\hspace{1em}(2)}$$
      由于$G_{n-1}=0$，则$A\beta =0$有有非零解，$\beta \ne 0$，对于第k行，

    $$({\varphi }_k,{\varphi }_0){\beta }_0+({\varphi }_k,{\varphi }_1){\beta }_1+...+({\varphi }_k,{\varphi }_{n-1}){\beta }_{n-1}=0$$

    即$({\varphi }_k,{\sum }_{i=0}^{n-1}{\beta }_i{\varphi }_i)=0,(k=0,1,...,n-1)$

    记$\Phi ={\sum }_{i=0}^{n-1}{\beta }_i{\varphi }_i$，我们可以得到
    $$(\Phi ,\Phi )=(\Phi ,\sum _{i=0}^{n-1}{\beta }_i{\varphi }_i)=\sum _{i=0}^{n-1}{\beta }_i(\Phi ,{\varphi }_i)=0$$
    根据函数内积的正定性，可以得到$\Phi =0$，即${\sum }_{i=0}^{n-1}{\beta }_i{\varphi }_i=0,\beta \ne 0$，得${\varphi }_j$线性相关，证毕。

    必要性一样可以用反证法证得！

    若函数系${\varphi }_i$线性相关，则存在不全为0的数$a_1,a_2,...,a_n$使得：

    $$\sum _{i=0}^{n-1}{a}_i{\varphi }_i=0$$

    等式两边同乘${\varphi }_1,{\varphi }_2,...,{\varphi }_n-1$，便得到方程组：$$\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_i({\varphi }_i,{\varphi }_j)=0$$

    因其有非零解，所以系数矩阵的行列式一定为0，与$G_{n-1}\ne 0$矛盾，得证。

函数的最佳平方逼近

• 定义：设 f ( x ) ∈ C [ a , b ] , Φ = s p a n { ϕ 0 , ϕ 1 , . . . , ϕ n } f(x) \in C[a,b], \Phi=span \{ \phi_0,\phi_1,...,\phi_n \} f(x)∈C[a,b],Φ=span{ϕ0​,ϕ1​,...,ϕn​}，若存在$S^{\ast }(x)\in \Phi$，使得

  $$||f-S^{\ast }|{|}_2^2=mi{n}_{S\in \Phi }||f-S|{|}_2^2=mi{n}_{S\in \Phi }{\int }_b^a\rho (x)[f(x)-S(x){]}^{2}dx$$

  称$S^{\ast }(x)$是$f(x)$在$\Phi \in C[a,b]$上的最佳逼近平方函数。

  其中，$S^{\ast }(x)$可以表示成$S^{\ast }(x)={\sum }_{j=0}^n{a}_j{\varphi }_j(x)$，等价于关于$(a_0,a_1,...,a_n)$函数的最小值：

  a r g m i n { I ( a 0 , a 1 , . . . , a n ) } = ∫ a b ρ ( x ) [ ∑ j = 0 n a j ϕ j ( x ) − f ( x ) ] 2 d x argmin\{ I(a_0,a_1,...,a_n) \} =\int_a^b \rho(x)[\sum_{j=0}^n a_j \phi_j(x)-f(x)]^2dx argmin{I(a0​,a1​,...,an​)}=∫ab​ρ(x)[j=0∑n​aj​ϕj​(x)−f(x)]2dx

  这是一个多元函数的极值问题，我们先利用其必要条件来试着求解，即在极值点处偏导数为0，

  $$\frac{\partial I}{\partial a_k}=2{\int }_a^b\rho (x)[\sum _{j=0}^n{a}_j{\varphi }_j(x)-f(x)]{\varphi }_k(x)dx=0,(k=0,1,...,n)$$

  可以得到称为法方程的约束方程：

  $$\sum _{j=0}^n{a}_j({\varphi }_k,{\varphi }_j)=(f,{\varphi }_k),(k=0,1,...,n)$$

  由于${\varphi }_j$线性无关，$G_n({\varphi }_0,{\varphi }_1,...,{\varphi }_n)\ne 0$，所以法方程有唯一解$a_k^{\ast }$，接下来证明所求解$S^{\ast }(x)=a_0^{\ast }{\varphi }_0(x)+a_1^{\ast }{\varphi }_1(x)+...+a_n^{\ast }{\varphi }_n(x)$为令$I$为最小值的解：

  即证明：${\int }_b^a\rho (x)[f(x)-S^{\ast }(x){]}^{2}dx\le {\int }_b^a\rho (x)[f(x)-S(x){]}^{2}dx$

  $$D={\int }_b^a\rho (x)[f(x)-S(x){]}^{2}dx-{\int }_b^a\rho (x)[f(x)-S^{\ast }(x){]}^{2}dx\text{\hspace{1em}(3)}$$

  这里我们需要用到$S^{\ast }(x)$为解的性质，即

  $$\sum _{j=0}^n{a}_j^{\ast }({\varphi }_k,{\varphi }_j)=(f,{\varphi }_k),(k=0,1,...,n)\to (f-\sum _{j=0}^n{\varphi }_j{a}_j^{\ast },{\varphi }_k)=0\to (f-S^{\ast }(x),{\varphi }_k)=0,(k=0,1,...,n)$$

  所以我们要凑出$(f-S^{\ast }(x),{\varphi }_k)$这一项，这一步稍显复杂，可以将式3化成下面的形式：

  $$D={\int }_b^a\rho (x)[S(x)-S^{\ast }(x){]}^{2}dx-2{\int }_b^a\rho (x)[S^{\ast }(x)-S(x)][f(x)-S^{\ast }(x)]dx\text{\hspace{1em}(4)}$$

  而$S^{\ast }(x)-S(x)$可以用$({\varphi }_0,{\varphi }_1,...,{\varphi }_n)$线性表出，所以式4中第二项积分为0，即有

  $$D={\int }_b^a\rho (x)[S(x)-S^{\ast }(x){]}^{2}dx\ge 0\text{\hspace{1em}(5)}$$

  即证明所求解$S^{\ast }(x)=a_0^{\ast }{\varphi }_0(x)+a_1^{\ast }{\varphi }_1(x)+...+a_n^{\ast }{\varphi }_n(x)$为令$I$为最小值的解

• 逼近误差

  令逼近误差$\delta =f-S^{\ast }(x)$，其平方误差为：

  $$||\delta |{|}_2^2=(f-S^{\ast },f-S^{\ast })=(f,f-S^{\ast })-(S^{\ast },f-S^{\ast })=(f,f)-(f,S^{\ast })=||f|{|}_2^2-\sum _{k=0}^n{a}_k^{\ast }(f,{\varphi }_k)\text{\hspace{1em}(6)}$$

  可以发现式6的第二项即法方程的右边项。

正交多项式

• 在求解最佳平方逼近时要法方程，为一个n阶线性方程组，由上可知，法方程存在唯一非零解，但是当n达到一定程度时，法方程的系数矩阵，即

  $$G_{n-1}=G({\varphi }_0(x),{\varphi }_1(x),...,{\varphi }_{n-1}(x))=\left|\begin{array}{cccc}({\varphi }_0,{\varphi }_0)& ({\varphi }_0,{\varphi }_1)& ...& ({\varphi }_0,{\varphi }_{n-1})\\ ({\varphi }_1,{\varphi }_0)& ({\varphi }_1,{\varphi }_1)& ...& ({\varphi }_1,{\varphi }_{n-1})\\ ⋮& ⋮& \ddots & ⋮\\ ({\varphi }_{n-1},{\varphi }_0)& ({\varphi }_{n-1},{\varphi }_1)& ...& ({\varphi }_{n-1},{\varphi }_{n-1})\end{array}\right|$$

  会是一个病态矩阵，求解结果受计算机舍入误差影响会剧烈抖动。

  这里如果当$G_n$为一个对角矩阵时，求解是十分简单直接的，不会存在抖动问题，这就要求${\varphi }_k$之间是互相正交的，即

  $$({\varphi }_i,{\varphi }_k)=0ifi!=kelse1$$

  同一线性空间的两组基是可以互相转换的，而将一组基转化为正交基有多种方法，最常见的便是施密特正交化，这在线代部分有介绍过给大家。

  其实很常见的一个多项式函数正交基就是$1,x,x^2,...,x^n$，利用施密特正交化给出的正交多项式$g_n(x)$有以下性质：

  1. $g_n(x)$是最高项系数为1的n次多项式
  2. 任一n次多项式$P_n(x)\in H_n$均可表示为$g_0(x),g_1(x),...,g_n(x)$的线性组合
  3. $g_n(x)$与任一次数小于n的多项式正交：很好理解，$g_n(x)$是$H_n$的一组基，则任何次数小于n的多项式都可以由$g_0,g_1,...,g_{n-1}$线性组成，则必与$g_n(x)$正交
  4. 存在$g_{n-1}(x)$、$g_n(x)$和$g_{n+1}(x)$的递推式：

  $$\begin{array}{rl}& g_{n+1}(x)=(x-{\alpha }_n)g_n(x)-{\beta }_n{g}_{n-1}(x),n=0,1,...,\\ & g_0(x)=1,g_{-1}(x)=0,\\ & {\alpha }_n=\frac{(x{g}_n,g_n)}{(g_n,g_n)}\\ & {\beta }_n=\frac{(g_n,g_n)}{(g_{n-1},g_{n-1})}\end{array}$$

  5. 设$g_0(x),g_1(x),...,g_n(x)$是在区间$[a,b]$上带权函数$\rho (x)$的正交多项式，则$g_n(x)$的n个根都是单重实根，且都在区间$[a,b]$

其他常见的正交多项式

Legendre（勒让德）多项式

• 区间为$[-1,1]$，权函数为1，其形式如下：

  P 0 ( x ) = 1 , P n ( x ) = 1 2 n n ! d n d x n { ( x 2 − 1 ) n } \large P_0(x)=1,P_n(x)=\frac{1}{2^n n!}\frac{d^n}{dx^n}\{ (x^2-1)^n \} P0​(x)=1,Pn​(x)=2nn!1​dxndn​{(x2−1)n}

  $${\int }_{-1}^1{P}_n(x)P_m(x)dx=\frac{2}{2n+1}ifn\ne melse0$$

  其中$(x^2-1{)}^n$是2n次多项式，则求导n次后，最高项系数为：

  $a_n=\frac{(2n)!}{2^n(n!{)}^2}$

  若将最高项系数归化为1，则归一化后的勒让德多项式形式为：

  P n ( x ) ‾ = n ! 2 n ! d n d x n { ( x 2 − 1 ) n } \large \overline{P_n(x)} =\frac{n!}{2n!}\frac{d^n}{dx^n}\{ (x^2-1)^n \} Pn​(x)​=2n!n!​dxndn​{(x2−1)n}

• 递推公式

  $(n+1)P_{n+1}(x)=(2n+1)x{P}_n(x)-n{P}_{n-1}(x),P_0(x)=1,P_1(x)=x$

Chebyshev（切比雪夫）多项式

• 区间为$[-1,1]$，权函数为$\rho (x)=\frac{1}{\sqrt{1-x^2}}$，其形式如下：

  $$T_n(x)=cos(n\arccos x),|x|\le 1$$

  如果令 $x=cos\theta$，则可以写成：

  $T_n(x)=\cos n\theta ,0\le \theta \le \pi$

• 递推公式：

  $$T_{n+1}(x)=2x{T}_n(x)-T_{n-1}(x),T_0(x)=1,T_1(x)=x$$

  由递推式我们可以知道，切比雪夫多项式的首项系数为$2^{n-1}$，且偶次项只有偶次幂，奇次项只有奇次幂。

函数按正交多项式展开

• Q: 若$f(x)\in C[a,b]$，用正交多项式 { g 0 , g 1 , . . . , g n } \{g_0,g_1,...,g_n\} {g0​,g1​,...,gn​}作基，求最佳平方逼近多项式

  $$S_n(x)=a_0{g}_0(x)+a_1{g}_1(x)+...+a_n{g}_n(x)$$

在前面最佳平方逼近内容中，我们得到法方程如下：

$$\sum _{j=0}^n{a}_j({\varphi }_k,{\varphi }_j)=(f,{\varphi }_k),(k=0,1,...,n)$$

由于 { g 0 , g 1 , . . . , g n } \{g_0,g_1,...,g_n\} {g0​,g1​,...,gn​}是正交多项式，可得
$$a_k=\frac{(f,g_k)}{(g_k,g_k)},(k=0,1,...,n)\text{\hspace{1em}(6)}$$

所以函数在正交多项式下的最佳平方逼近为：

$$f\to S_n(x)=\sum _{k=0}^n\frac{(f,g_k)}{(g_k,g_k)}{g}_k(x)$$

注意使用勒让德或者切比雪夫多项式逼近函数时要在定义域$[-1,1]$内，如果不在此范围，需要用线性变换使x的范围落在$[-1,1]$内。

若$x\in [a,b]$，则令$x=\frac{b-a}{2}t+\frac{b+a}{2},-1\le t\le 1$，即$F(t)=f(\frac{b-a}{2}t+\frac{b+a}{2})$，求得对应的最佳平方逼近函数$S^{\ast }(x)$后，$f(x)$对应的最佳平方逼近函数即为$S^{\ast }(\frac{1}{b-a}(2x-a-b))$

正交函数作最小二乘拟合

• 给出对应的点集 { x i } , { y i } \{ x_i \}, \{y_i \} {xi​},{yi​}，取带权函数$\omega (x)$的正交函数集为 { ϕ 0 , ϕ 1 , . . . , ϕ n } \{ \phi_0,\phi_1,...,\phi_n \} {ϕ0​,ϕ1​,...,ϕn​}，注意这里的点是离散的，而我们到上面为止所提到的内积都是基于函数定积分，离散的正交拟合在计算机上实现上更便捷。

• 由式6可知，离散形式的各正交基系数为：

  $$a_k^{\ast }=\frac{(f,{\varphi }_k)}{({\varphi }_k,{\varphi }_k)}=\frac{{\sum }_{i=0}^m\omega (x_i)f(x_i){\varphi }_k(x_i)}{{\sum }_{i=0}^m\omega (x_i){\varphi }_k^2(x_i)},(k=0,1,...,n)$$

  其实正交多项式最小二乘拟合与常见的最小二乘拟合（在之前矩阵论中我们曾介绍过）很类似，但要注意区间范围和权函数等。理论上正交多项式所构成的系数矩阵避免了病态问题，求解是更加高效的。

代码实现

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches


class Orthogonal_Polynomials_Fitting:
    """ Fitting points by Orthogonal polynomials

    Author: Junno
    Date: 2022-03-24
    """

    def __init__(self, x, y):
        """__init__ 

        Args:
            x ([np.array]): [x of points]
            y ([np.array]): [y of points]
        """
        assert x.shape[0] == y.shape[0] or len(x) == len(y)
        self.x = x
        self.y = y
        self.N = self.x.shape[0]
        self.LP_records = {}
        self.CH_records = {}
        self.low = np.min(x)
        self.high = np.max(x)
        self.shift = False
        # if x's range not in [-1,1], do linear shift
        if self.low < -1 or self.high > 1:
            self.shift = True
            self.shift_x = np.arange(-1, 1, 2./self.N)

    @classmethod
    # recursion formula of Legender Polynomials
    def Legender_Polynomials(self, x, n):
        if n == 0:
            return 1
        elif n == 1:
            return x
        else:
            if (x, n) in self.LP_records.keys():
                return self.LP_records[(x, n)]
            else:
                ans = ((2*n-1)*x*Legender_Polynomials(x, n-1) -
                       (n-1)*Legender_Polynomials(x, n-2))/(n)
                self.LP_records[(x, n)] = ans
                return ans

    @classmethod
    # recursion formula of Chebyshev Polynomials
    def Chebyshev_Polynomials(self, x, n):
        if n == 0:
            return 1
        else:
            # let x=cos(theta), 0<=theta<=pi, Tn(x)=cos(n*theta)
            return np.cos(n*x)
#         elif n==1:
#             return x
#         else:
#             if (x,n) in self.CH_records.keys():
#                 return self.CH_records[(x,n)]
#             else:
#                 ans=2*x*Chebyshev_Polynomials(x,n-1)-Chebyshev_Polynomials(x,n-2)
#                 self.CH_records[(x,n)]=ans
#                 return ans

    def calculate(self, order, method='Legender', show_error=True, show_fitting_result=True):
        """calculate_I calculate y for order and method

        Args:
            order ([int], optional): [interpolation order]. Defaults to 1.
            show_error (bool, optional): [print fitting error]. Defaults to True.
            show_fitting_result (bool, optional): [show fitting result by plt]. Defaults to True.

        Returns:
            [predict y]: [as show]
        """
        cx = self.shift_x if self.shift else self.x

        if method == 'Legender':
            Method = Legender_Polynomials
        elif method == 'Chebyshev':
            Method = Chebyshev_Polynomials
            # with x=cosθ,we can let the weight function to be 1
            cx = np.arange(0, np.pi, np.pi/cx.shape[0])

        A = np.zeros((order, order))
        for i in range(order):
            for x_ in cx:
                A[i, i] += Method(x_, i)*Method(x_, i)

        B = np.zeros((order, 1))
        for i in range(order):
            for x_, y_ in zip(cx, self.y):
                B[i, 0] += Method(x_, i)*y_

        # solute coefficients
        ai = np.linalg.inv(A)@B

        predict_y = np.zeros_like(self.y)
        for j in range(order):
            for i in range(cx.shape[0]):
                #                 x_ = self.x[i] if not self.shift else 1/(self.high-self.low)*(2*self.x[i]-self.high-self.low)
                predict_y[i] += Method(cx[i], j)*ai[j, 0]

        error=np.linalg.norm(y-predict_y)
        
        if show_error:
            print("Fitting Error: {:.5f}".format(error))

        if show_fitting_result:
            plt.subplot(111)
            plt.scatter(self.x, self.y, s=5, marker='*', c='green')
            plt.scatter(self.x, predict_y, s=5, marker='o', c='red')
            Raw = mpatches.Patch(color='green', label='Raw Points')
            Predict = mpatches.Patch(color='red', label='Predict Points')
            plt.title('Orthogonal-Polynomials by {} with {} orders'.format(method,order))
            plt.legend(handles=[Raw, Predict], loc='best')
            plt.show()
        return predict_y,error

Test example

# prepare some point set
F=lambda x: np.cos(2*np.exp(x))
x=np.arange(1,2,0.02)
y=F(x)
OP=Orthogonal_Polynomials_Fitting(x,y)

# Legender Method
py, _=OP.calculate(5,method='Legender')
>>> Fitting Error: 1.50786

py, _=OP.calculate(6,method='Legender')
>>> Fitting Error: 0.91998

# Chebyshev Method
py, _=OP.calculate(5,method='Chebyshev')
>>> Fitting Error: 1.45608

py, _=OP.calculate(6,method='Chebyshev')
>>> Fitting Error: 1.02585

# show error with order and methods
LG_Error_with_order=[]
CH_Error_with_order=[]
for i in range(1,11):
    py,error=OP.calculate(i,method='Legender',show_error=False,show_fitting_result=False)
    LG_Error_with_order.append(error)
    py,error=OP.calculate(i,method='Chebyshev',show_error=False,show_fitting_result=False)
    CH_Error_with_order.append(error)

plt.subplot(111)
plt.plot(np.arange(0,10,1),LG_Error_with_order, marker='*',c='green')
plt.plot(np.arange(0,10,1),CH_Error_with_order,marker='o',c='red')
LG_Error=mpatches.Patch(color='green', label='LG_Error_with_order')
CH_Error=mpatches.Patch(color='red', label='CH_Error_with_order')
plt.title('Errors with different Orthogonal-Polynomials and Orders')
plt.legend(handles=[LG_Error,CH_Error],loc='best')
plt.show()

最佳一致逼近

基础内容

• 前面我们提到一致逼近或者说均匀逼近的误差衡量形式为：$||f(x)-P(x)|{|}_{\infty }=ma{x}_{a<x<b}|f(x)-P(x)|$，一般我们做函数逼近时比较喜欢用均方逼近，即最小二乘的误差衡量，因为平方的引入使得误差函数容易求导，而最佳一致逼近则是一种基于绝对值的最小最大问题，不易求导，研究起来稍微会困难一些。但是，没有跨不过的困难，只有勇敢的Chebyshev！

• 最佳一致逼近多项式OUAP（optimal uniform approximating polynomial）

  设 f ∈ C [ a , b ] , p n ∈ H n = s p a n { 1 , x , . . . , x n } f \in C[a,b], p_n \in H_n=span\{ 1,x,...,x^n \} f∈C[a,b],pn​∈Hn​=span{1,x,...,xn}，记$\Delta (f,p_n)=||f-p_n|{|}_{\infty }=ma{x}_{a\le x\le b}|f(x)-p_n(x)|$为$p_n$与$f$的偏差。由定义知，$\Delta (f,p_n)\ge 0$。

  若$x_0$满足$p(x_0)-f(x_0)=\pm ||P-f|{|}_{\infty }$，称$x_0$为正or负偏差点。

  在$H_n$中，如果存在满足，$E_n=in{f}_{p_n\in H_n}\Delta (f,p_n)=mi{n}_{p_n\in H_n}\Delta (f,p_n)=mi{n}_{p_n\in H_n}[ma{x}_{a\le x\le b}|f(x)-p_n(x)|]$，称$E_n$为$p_n$与$f$的最小偏差，满足上式的$p_n$称为$f(x)$在$[a,b]$上的最佳一致逼近多项式。

• 若$P(x)\in H_n$是OUAP，则同时存在正负偏差点

  证明如下：采用反证法，设只有正偏差点

  对于所有$x\in [a,b]$，都有$P(x)-f(x)>-E_n$，因$P(x)-f(x)$在$[a,b]$连续，故最小值大于$-E_n$，记为$-E_n+2h,h>0$,即有：

  $-E_n+2h\le P(x)-f(x)\le E_n$

  $-E_n+h\le [P(x)-h]-f(x)\le E_n-h$

  $|[P(x)-h]-f(x)|\le E_n-h$

  则可得$P(x)-h$与$f(x)$的偏差小于$E_n$，与给定条件矛盾，得证OUAP同时存在正负偏差点。

• OUAP的充分必要条件（Chebyshev 最佳一致逼近定理）

  $P(x)\in H_n$是$f(x)\in C[a,b]$的最佳一致逼近多项式的充要条件为$P(x)$在$[a,b]$上至少有n+2个轮流的正负偏差点，即有n+2个点$a\le x_1<x_2<...<x_{n+2}\le b$，使得
  $$P(x_k)-f(x_k)=(-1{)}^k\sigma ||P(x)-f(x)|{|}_{\infty },\sigma =\pm 1,k=1,2,...,n+2$$

  图示如下：

  

• OUAP的存在性与唯一性

  • 存在性证明：Borel定理，证明复杂从略

  • 唯一性证明：

    采用反证法证明，设存在两个最佳一致逼近多项式$P(x)$和$Q(x)$，则对于$x\in [a,b]$，都有：

    $-E_n\le P(x)-f(x)\le E_n$，$-E_n\le Q(x)-f(x)\le E_n$

    可以简单验证得到对于P和Q的加权和为1的线性组合都是$H_n$的OUAP。

    将上面两个式子各以$1/2$加权然后相加，得：

    $-E_n\le \frac{P(x)+Q(x)}{2}-f(x)\le E_n$

    则$R(x)=\frac{P(x)+Q(x)}{2}$也是OUAP，有$n+2$个交错点 { x k } \{ x_k \} {xk​}满足$R(x_k)-f(x_k)=(-1{)}^k\sigma E_n,(k=1,2,...,n+2)$

    由偏差点的定义，可得：

    $$E_n=|R(x_k)-f(x_k)|=|\frac{P_{(}{x}_k)-f(x_k)}{2}+\frac{Q_{(}{x}_k)-f(x_k)}{2}|\text{\hspace{1em}(7)}$$

    由于$|P(x_k)-f(x_k)|\le E_n,|Q(x_k)-f(x_k)|\le E_n$，当且仅当

    $$\frac{P_{(}{x}_k)-f(x_k)}{2}=\frac{Q_{(}{x}_k)-f(x_k)}{2}=\pm \frac{E_n}{2}$$

    时式7成立，即有$P_{(}{x}_k)-f(x_k)=Q_{(}{x}_k)-f(x_k)$，得$P(x_k)=Q(x_k)$。说明$P(x)-Q(x)\in H_n$有n+2个零点，只能得到$P(x)-Q(x)=0$，即$P(x)=Q(x)$，与条件设定矛盾，证得OUAP是唯一的。

如何得到OUAP

• 从上面的内容我们可以得到，OUAP有n+2个正负交错点，即有n+1个根，所以$P_n(x)$可以看成是$f(x)$的一个n阶插值多项式。 以Lagrange插值法为例，那我们接下来的问题就是确定这n+1个插值节点的位置，使得得到的插值函数正好是$f(x)$的OUAP。

• 根据OUAP的条件和Lagrange的误差余项，满足OUAP的插值多项式即为其对应的误差余项的无穷范数最小，即

  $$||R_n(x)|{|}_{\infty }=||\frac{f^{n+1}(\xi )}{(n+1)!}\prod _{i=0}^n(x-x_i)|{|}_{\infty }$$

  其中，假设$f^{n+1}(\xi )$连续且有上限，则我们要求 { x 1 , . . , x n } \{x_1, .., x_n \} {x1​,..,xn​}使得$w_n(x)={\prod }_{i=0}^n(x-x_i)$形式下的$||w_n(x)|{|}_{\infty }$最小。（注意，这里我们的约束条件是$w_n(x)$的上界最小，一不定就是最小的$w(x)$，所以这是一个近似求解）

  假设$w_n(x)$是最高项系数为1的n阶多项式，可以表示为$w_n(x)=x^n-P_{n-1}(x)$，则问题转变为求$x^n$的n-1阶OUAP。

  这里求得$x^n$的n-1阶OUAP不是通过直接求解，前面定理说明OUAP是存在且唯一的，即只要我们找到一个满足条件的OUAP（存在n+1个交错正负偏差点），那就可以解决问题。

• （虽然感觉很突兀，但老师也是这么说的）考虑$[-1,1]$上Chebyshev多项式，$T_n(x)=\cos (n\arccos (x))$，在$x_k=\cos (\frac{k}{n}\pi ),(k=0,1,...,n)$时，$T_n(x)$交错取得极大值1和极小值-1；当$x_k=\cos (\frac{2k-1}{2n}\pi ),(k=1,...,n)$时，$T_n(x)=0$，即有n个零点。

  可以令归一化的Chebyshev多项式（最高项系数为1）为$\frac{1}{2^{n-1}}{T}_n(x)=x^n-P_{n-1}^{\ast }(x)$，即其无穷范数可以表示为
  $$||\frac{1}{2^{n-1}}{T}_n(x)|{|}_{\infty }=\frac{1}{2^{n-1}}||T_n(x)|{|}_{\infty }=\frac{1}{2^{n-1}}$$

• 结论

  即在$[-1,1]$上求使得$w_n(x)={\prod }_{i=0}^n(x-x_i)$的无穷范数最小的形式即为归一化的Chebyshev多项式，要求的插值节点即为$T_n(x)$的n+1个零点。

  对于一般区间的$x\in [a,b]$，可作变量替换$x=\frac{a+b}{2}+\frac{b-1}{2}t$，则$t\in [-1,1]$，此时有：

  $$w_{n+1}(x)=w_{n+1}(\frac{a+b}{2}+\frac{b-a}{2}t)=(\frac{a+b}{2}+\frac{b-1}{2}t-x_0)...(\frac{a+b}{2}+\frac{b-1}{2}t-x_n)=(\frac{b-a}{2}{)}^{n+1}(t-t_0)...(t-t_n)=(\frac{b-a}{2}{)}^{n+1}\frac{1}{2^n}{T}_{n+1}(t)$$

  即以$x_k=\frac{a+b}{2}+\frac{b-a}{2}\cos (\frac{2k+1}{2n+2})\pi$为插值节点，$(k=0,1,...,n)$。其余项：

  $R_n(x)=\frac{f^{n+1}(\xi )}{(n+1)!}\frac{(b-a{)}^{n+1}}{2^{2n+1}}{T}_{n+1}(t)$有最小上界。

代码实现

• 以$f(x)=e^x$为例， 若控制误差余项在$0.5\times 10^{-4}$内，则Chebyshev插值节点的个数n可由$|R_n(x)|\le \frac{e}{(n+1)!}\times \frac{1}{2^{2n+1}}<\frac{1}{2}\times 10^{-4}$确定，即n=4。下面验证以上过程及最后的误差余项上界是否满足条件：

# https://blog.csdn.net/weixin_40519529/article/details/123215263
import numpy as np
from functools import reduce
import matplotlib.pyplot as plt

# Lagrange 线性插值函数
class Lagrange_Polynomial_Interpolation(object):
    """Lagrange_Polynomial_Interpolation 

    base function: Ln(x)=\sum_{k=0}^{n} yk l_k(x)
    Author: Junno
    Date: 2022-03-01
    """

    def __init__(self, x, y):
        """__init__ 

        Args:
            x ([np.array]): [x of points]
            y ([np.array]): [y of points]
        """
        self.N = x.shape[0]
        self.x = np.copy(x)
        self.y = np.copy(y)

    def calculate_lk(self, x, k,):
        """calculate_lk(x)

        base function: lk(x)=∏——{j=0,j!=k}^{n} (x-xj)/(xk-xj)
        Args:
            x ([float]): [x]
            k ([int]): [k]

        Returns:
            [float]: [lk(x)]
        """
        ta = reduce(lambda x, y: x*y, [x-self.x[i]
                    for i in range(self.N) if i != k])
        tb = reduce(lambda x, y: x*y, [self.x[k]-self.x[i]
                    for i in range(self.N) if i != k])
        return ta/tb

    def predict(self, x):
        """predict 

        Args:
            x ([float]): [x]

        Returns:
            [float]: [predict y]
        """
        return np.sum([self.y[k]*self.calculate_lk(x, k) for k in range(self.N)])

# points
x0=0
xn=1
x=np.arange(x0,xn,1/10)
F=lambda x :np.exp(x)
y=F(x)

order=5
Cheby_zero_points=np.array([np.cos(np.pi/(2*order)*(2*i+1)) for i in range(order)])
nx=(x0+xn)/2+(xn-x0)/2*Cheby_zero_points
nx=np.sort(nx)
print(nx)
>>> array([0.02447174, 0.20610737, 0.5       , 0.79389263, 0.97552826])

LP=Lagrange_Polynomial_Interpolation(nx,F(nx))
py=np.array([LP.predict(k) for k in x])
print(np.max(np.abs(py-y)))
>>> 2.885868676028025e-05

• 如果使用不同数量的Chebyshev插值节点来构造，最后得到的误差上限变化如下图：

  

未完待续

• 码字不易，给个赞吧~