数值分析1--拉格朗日插值法

原理介绍

1、定义

2、满足插值条件的Ln(x)是否唯一？
假设Ln(x)为n次多项式，则带入各插值点条件可以得到如下n+1个方程：

我们要证Ln(x)的唯一性，就是证a0到an的唯一性。可以将1，x0，x0^2…当作系数矩阵，其系数行列式如下：

很明显这是一个范德蒙行列式，根据其计算方法，知该行列式的解为：

又因为x选自坐标系曲线上不同的点，因此若i != j，则x_i != x_j,则该行列式不为0。
最后由非齐次方程组AX = B的解唯一的条件是R(A)=R(A,B)=n,很明显满足，因此：

3、待定系数法求Ln(x)
如下图所示，唯一需要注意的是题目给出了三个零点，两个导数零点，因此一共是五个条件，我们知道五个条件可以解出五元函数，即a0-a4，此时我们所设的多项式应该是个四次多项式，其它内容就是小学数学。

4.拉格朗日多项式
定义：
描述：将n+1个y_i作为多项式的各个系数。

举例：写出了L1(x)的点斜式，进而转换成两点式，成功转换成上述形式。

n >= 1时：

例：

5.拉格朗日多项式插值余项

定义为插值余项，也就是截断误差。
用罗尔定理及其推广得到下面的定义式。

6.拉格朗日通项公式如下：

7.quiz

给了六个插值点，则根据l(x)的通项表达式可以知道l(x)是一个五次多项式曲线，如A所示。

代码实现

拉格朗日类的定义：

import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
from Lagrangeinterpolation.utils import interp_utils


class LagrangeInterpolation:
    def __init__(self, x, y):
        self.x = np.asarray(x, dtype=float)  # as_array是视图， array是副本
        self.y = np.asarray(y, dtype=float)
        if len(self.x) > 1 and len(self.x) == len(self.y):  # 健壮性判断，要知道自己要判断什么，数据传输不能有误
            self.n = len(self.x)
        else:
            raise ValueError("插值数据(x, y)维度不匹配")
        self.polynomial = None  # 最终多项式
        self.poly_coefficient = None  # 多项式系数
        self.coefficient_order = None  # 多项式系数对应地次数
        self.y0 = None  # 所求插值点的值

    def fit_interp(self):
        """
        核心算法：生成拉格朗日多项式
        :return:
        """
        x = sym.Symbol('x')
        self.polynomial = 0.0  # 插值多项式实例化
        for i in range(self.n):  # 传入10个插值点，0-9，len（x）就是10，所以遍历所有插值点就是range（self.n）
            basic_fuc = self.y[i]
            for j in range(i):
                basic_fuc *= (x - self.x[j]) / (self.x[i] - self.x[j])
            for j in range(i + 1, self.n):
                basic_fuc *= (x - self.x[j]) / (self.x[i] - self.x[j])
            self.polynomial += basic_fuc
        self.polynomial = sym.expand(self.polynomial)
        polynomial = sym.Poly(self.polynomial, x)
        self.poly_coefficient = polynomial.coeffs()
        self.coefficient_order = polynomial.monoms()

    def cal_interp_x0(self, x0):
        self.y0 = interp_utils.cal_interp_x0(self.polynomial, x0)
        return self.y0

    def plt_interpolation(self, x0, y0):
        params = (self.polynomial, self.x, self.y, 'Lagrange', x0, y0)
        interp_utils.plt_interpolation(params)

工具函数：

import numpy as np
import matplotlib.pyplot as plt


def cal_interp_x0(polynomial, x0):
    """
    计算给定插值点的数值，即插值
    :param polynomial:
    :param x0:插值横坐标
    :return:
    """
    x0 = np.asarray(x0, dtype=np.float64)
    n0 = len(x0)  # 所求插值点的个数
    y_0 = np.zeros(n0)  # 存储插值点对应的插值
    t = polynomial.free_symbols.pop()
    for i in range(n0):
        y_0[i] = polynomial.evalf(subs={t: x0[i]})
    return y_0


def plt_interpolation(params):
    polynomial, x, y, title, x0, y0 = params
    plt.figure(figsize=(8, 6))
    plt.plot(x, y, "yo", label="Interpolation base point")
    xi = np.linspace(min(x), max(x), 100)
    yi = cal_interp_x0(polynomial, xi)
    plt.plot(xi, yi, "r--", label="Interpolation polynomial")
    if x0 is not None and y0 is not None:
        plt.plot(x0, y0, "g*", label="Interpolation point values")
    plt.legend()
    plt.xlabel("x", fontdict={"fontsize": 12})
    plt.ylabel("y", fontdict={"fontsize": 12})
    plt.title(title + "interpolation polynomial and values")
    plt.grid(ls=":")
    plt.show()

主函数

from Lagrangeinterpolation.lagrange import LagrangeInterpolation
import numpy as np

x = np.linspace(0, 2 * np.pi, 10, endpoint=True)
y = np.sin(x)
x0 = np.array([1.2, 1.3, 1.5, 2, 3, 4, 5, 6])
lag_interp = LagrangeInterpolation(x=x, y=y)
lag_interp.fit_interp()
print("拉格朗日多项式如下：")
print(lag_interp.polynomial)
print("拉格朗日插值多项式系数向量和对应阶次：")
print(lag_interp.poly_coefficient)
print(lag_interp.coefficient_order)
y0 = lag_interp.cal_interp_x0(x0)
lag_interp.plt_interpolation(x0, y0)

结果展示

摘录于：（1）https://www.bilibili.com/video/BV1Lq4y1U7Hj/?spm_id_from=333.337.search-card.all.click
（2）https://www.bilibili.com/video/BV1AK4y1k7Px/?spm_id_from=333.337.search-card.all.click&vd_source=35775f5151031f11ee2a799855c8e368