拉格朗日插值


Ⅰ. 作用 –> Ⅱ. 原理 –> Ⅲ. 余项 –> Ⅳ. 示例 –> Ⅴ. Python 代码 –> Ⅵ. 背景


Ⅰ. 作用

顾名思义，既然是插值，那就可以补充数据，在旧数据之间插入新值. 拉格朗日插值得到的多项式当然也可用于预测，但我强烈不建议你去预测那些遥远的节点，一来，变数太大不够准确，二来，插值构造的多项式阶数过高，他们说，这会产生龙格现象.


Ⅱ. 原理

已知从小到大排序的唯一点列，统计上亦称“观测”，observation：

\overrightarrow{X}=[x_1,x_2,\dots,x_n]^{\top}

\overrightarrow{Y}=[y_1,y_2,\dots,y_n]^{\top}

定义拉格朗日插值基函数：

\forall t\in\{1,2,\dots,n\},\ \ \ \ l_{t}(x)=\dfrac{\prod\limits_{i=1}^{n\setminus t}(x-x_{i})}{\prod\limits_{i=1}^{n\setminus t}(x_t-x_{i})}

定义拉格朗日插值多项式：

\begin{gather*} L(x)=\sum\limits_{t=1}^{n}y_tl_{t}(x) \end{gather*}

阶数：

\partial l_{t}(x)=n-1,\ \ \ \ \partial L(x)=n-1

新的插值节点：

x^{*}\Longrightarrow y^{*}=L(x^{*})


Ⅲ. 余项

插值多项式毕竟不是本源函数（一般现实里本源函数也是未知），因此，拉格朗日插值多项式就存在误差，亦称余项.

假定人类已知性质良好的本源函数：

y=f(x)

则，拉格朗日插值多项式余项（截断误差）为：

\{x_1,x_2,\dots,x_n\}\subseteq[a,b]\ \ \ \&\ \ \ \forall k,\ x_k<x_{k+1},\ y_k=f(x_k)

{\rm Residual}(x)=f(x)-L(x)=\dfrac{f^{n}(\xi)}{n!}\times\prod\limits_{i=1}^{n}(x-x_i),\ \ \ \xi\in(a,b)


Ⅳ. 示例

已知观测：

\overrightarrow{X}=[1,2,3]^{\top}

\overrightarrow{Y}=[1,-4,9]^{\top}

则二次拉格朗日插值多项式为：

L(x) &= 1\times\dfrac{(x-2)(x-3)}{(1-2)(1-3)}-4\times\dfrac{(x-1)(x-3)}{(2-1)(2-3)} \\ &+ 9\times\dfrac{(x-1)(x-2)}{(3-1)(3-2)}

插值：

x^{*}=1.5

y^{*}=L(x^{*})=-3.75

预测遥远的节点：

x^{*}=12

y^{*}=L(x^{*})=936

注意到，由已知点列构造的拉格朗日插值多项式，对于新的插值节点 12，距离自变量 1, 2, 3 整体较远，代入得到的函数值 936 非常大. 用插值多项式预测遥远的节点，变数极大，结果是不稳定的.


Ⅴ. 代码

Lagrange.py
  1'''
  2# System --> Windows & Python3.8.0
  3# File ----> Lagrange.py
  4# Author --> Illusionna
  5# Create --> 2023/11/02 15:53:15
  6'''
  7# -*- Encoding: UTF-8 -*-
  8
  9
 10import copy
 11import numpy as np
 12from typing import Literal
 13from sympy import (symbols, expand)
 14
 15
 16class LAGRANGE:
 17    """
 18    --------------------------------------------------
 19    类: 拉格朗日插值多项式.
 20        obj = LAGRANGE(X:list, Y:list)
 21    --------------------------------------------------
 22    函数:
 23    1. obj.Interpolate(x:float) -> float
 24        x 为插值节点.
 25    2. obj.Show(precision:int=12, mode='definition') -> None
 26        precision 为显示定义式的多项式系数的精度;
 27        mode 为多项式显示的模式 --> 'simplify' & 'definition'
 28    --------------------------------------------------
 29    """
 30    def __init__(
 31        self,
 32        *args,
 33        X:list,
 34        Y:list,
 35        **kwargs
 36    ) -> None:
 37        self.X = X
 38        self.Y = Y
 39        self.__base = LAGRANGE.__BaseCoefficients(self)
 40
 41    def Interpolate(self, x:float) -> float:
 42        """
 43        拉格朗日插值.
 44        """
 45        result = 0
 46        val = x
 47        for i in range(0, len(self.X), 1):
 48            temp = list(
 49                map(lambda x: val - x, self.X)
 50            )
 51            temp.pop(i)
 52            numerator = np.array(temp).prod()
 53            del temp
 54            # --------------------------------------
 55            """
 56            如果想获取更精确的插值，解锁如下注释...
 57            """
 58            # temp = list(
 59            #     map(lambda x: self.X[i] - x, self.X)
 60            # )
 61            # temp.remove(0)
 62            # denominator = np.array(temp).prod()
 63            # del temp
 64            """
 65            用如下注释顶替 result 输出结果...
 66            """
 67            # result = result + (self.Y[i] * numerator / denominator)
 68            # --------------------------------------
 69            result = result + self.__base[i]*numerator
 70        return result
 71
 72    def Show(
 73        self,
 74        precision:int=12,
 75        mode:Literal['definition', 'simplify']='definition'
 76    ) -> None:
 77        """
 78        控制台显示拉格朗日多项式.
 79        """
 80        if mode == 'definition':
 81            showString = '\033[036mL(x)\033[0m = '
 82            for i in range(0, len(self.__base), 1):
 83                coef = self.__base[i]
 84                string = LAGRANGE.__PolynomialString(self.X, i, 'definition')
 85                temp = f'\033[033m%.{precision}f\033[0m{string} \033[031m+\033[0m ' % coef
 86                showString = showString + temp
 87            showString = showString[:-13]
 88            del temp
 89            print(showString)
 90        elif mode == 'simplify':
 91            showString = ''
 92            for i in range(0, len(self.__base), 1):
 93                coef = self.__base[i]
 94                string = LAGRANGE.__PolynomialString(self.X, i, 'simplify')
 95                string = string[:-1]
 96                temp = f'%.{precision}f*{string}+' % coef
 97                showString = showString + temp
 98            showString = showString[:-1]
 99            temp = str(expand(showString))
100            expression= 'L(x) = '
101            expression = expression + temp
102            del temp
103            print(expression)
104        else:
105            print('Error...')
106            exit(0)
107
108    def __BaseCoefficients(self) -> list:
109        coefficientsVector = []
110        for i in range(0, len(self.Y), 1):
111            y = self.Y[i]
112            temp = list(
113                map(lambda x: self.X[i] - x, self.X)
114            )
115            temp.remove(0)
116            denominator = np.array(temp).prod()
117            coefficientsVector.append(y / denominator)
118        del temp
119        return coefficientsVector
120
121    def __PolynomialString(vector:list, i:int, mode:str) -> str:
122        temp = copy.deepcopy(vector)
123        temp.pop(i)
124        string = ''
125        if mode == 'definition':
126            for j in range(0, len(temp), 1):
127                value = temp[j]
128                if value > 0:
129                    string = string + f'(x-{value})'
130                elif value < 0:
131                    string = string + f'(x+{abs(value)})'
132                elif value == 0:
133                    string = string + '(x)'
134            del temp
135            return string
136        else:
137            for j in range(0, len(temp), 1):
138                value = temp[j]
139                if value > 0:
140                    string = string + f'(x-{value})*'
141                elif value < 0:
142                    string = string + f'(x+{abs(value)})*'
143                elif value == 0:
144                    string = string + '(x-0)*'
145            del temp
146            return string
147
148
149if __name__ == '__main__':
150    """
151    以 y = (x^4)*(e^x) 为例.
152    查看 LAGRANGE 类文档
153    >>> print(LAGRANGE.__doc__)
154    """
155    # 测试拉格朗日插值类效果.
156    print('\033[H\033[J', end='')
157    print(LAGRANGE.__doc__)
158
159    # ----------------------------------------
160    # 插值核心代码.
161    X = [-7, -6.2, -5.4, -4.6, -3.8, -3]
162    Y = [2.18, 2.99, 3.84, 4.50, 4.66, 4.03]
163    obj = LAGRANGE(X=X, Y=Y)
164    value = obj.Interpolate(-5)
165    # ----------------------------------------
166
167    print(f'当 x = -5, 插值 L(x) = {value}')
168    print('\n插值结果定义式:')
169    obj.Show(precision=7, mode='definition')
170    print('\n插值结果化简式')
171    obj.Show(mode='simplify')
172    print('')
173
174    # ----------------------------------------
175
176    import matplotlib.pyplot as plt
177
178    x = np.linspace(-7, -1, 20)
179    y1 = x**4 * np.exp(x)
180    y2 = []
181    for i in range(0, len(x), 1):
182        y2.append(obj.Interpolate(x[i]))
183
184    observation = plt.plot(X, Y, 'bo')
185    interpolation = plt.plot(x, y2, 'r*')
186    function = plt.plot(x, y1, 'g-')
187
188    plt.title('Lagrange Interpolation')
189    plt.legend(['observation', 'interpolation', 'function: $y=x^4e^x$'])
190    plt.show()

插值结果：

figure

Ⅵ. 背景

【背景根据历史虚构】

在 16 世纪，由于欧洲航海事业的发展，先驱们就必须掌握“观天象”这一本领，自然就离不开天文学的引导.

十七世纪新纪元的伊始之年，步 Giordano Bruno 火刑之后尘，欧洲的人们逐渐接受《日心说》思想.

假设在 17 世纪，人类通过某些特殊手段获取到太阳和地球的部分数据，包括距离和引力，但由于马虎，忘了记载数量级、量纲单位 :)>

先驱们留下的资料

距离（Distance）

引力（Gravity）

1.46

3.7160276576

1.47

3.6656414249

1.48

3.6162730803

1.50

3.5204820244

1.51

3.4740075238

1.52

3.4284472624

人们知道一年有 365 天、四个季度，地球围绕太阳椭圆公转，因此，地日距离自然不同，表格中六个观测的距离是地表最强肉眼观察大师第谷所得，（六个引力值是未来的笔者我所提供）.

某一天，第谷又观察到一个新的地日距离 1.49 无单位，人们想知道当天的引力值，但未来的我因网络信号不佳，提供的引力值传输到数据链路层被吃掉了.

焦头烂额之际，主人公 Joseph-Louis Lagrange 粉墨登场，说：“我会”，于是他创造性地给出拉格朗日插值法，熟练 Lagrange.py 代码，在 1.48 和 1.50 之间插入 1.49 的函数值，终端返回：

3.5678953899295365

至此，人们的问题得以解决.

但不久之后，质疑的声音在人群中嘹响. 譬如，未来的我，我觉得拉格朗日的结果不一定对，他方法得到的结果应该有很大偏差. 所以，这个偏差究竟该怎么断定？

这个时候，我们就需要这六个观察对应的本源函数，用绝对真值和估计插值比对，看看误差大小是否在我们可控范围之内.

纵使世间万象万千，好在稍后的牛顿发现了现象背后的本源函数，即万有引力定律：

\begin{gather*} F=\dfrac{GM_1M_2}{r^2}\Longrightarrow {\rm Gravity(Distance)}=\dfrac{GM_1M_2}{{\rm Distance}^2}=\dfrac{k}{r^2} \end{gather*}

通过表格已知的数据，把分子作为整体消去后，我们得到距离为 1.49 时的万有引力：

3.5678953898473048

通过本源函数计算的引力值和拉格朗日插值得到的引力值，一直到小数点后 10 位才发生分歧（不同电子设备计算的值可能略有不同，此处只表意），面对如此穷极至小的误差，我拜倒在拉格朗日的膝底.

从此，笔者我质疑的声音消失了.

当然，消去分子整体后，重新确定万有引力函数的常数项系数，上述过程，通过拉格朗日插值多项式余项估算也可以的：

{\rm Residual}(1.49)=F(1.49)-L(1.49)=\dfrac{\dfrac{\partial^6 F}{\partial r}\bigg|_{r=\xi}}{720}\times\prod\limits_{i=1}^{6}(1.49-r_i),\ \ \ \xi\in(1.46,1.52)