9.2 拉格朗日插值

醒过来摸鱼

已于 2022-06-28 23:25:19 修改

阅读量1.3k

点赞数

分类专栏：数学算法

于 2022-04-04 11:40:41 首次发布

本文链接：https://blog.csdn.net/m0_66201040/article/details/123950255

版权

拉格朗日插值多项式 Python实现化简霍纳法则

关键词由CSDN通过智能技术生成

数学算法专栏收录该内容

42 篇文章 16 订阅

订阅专栏

拉格朗日插值的公式是
$P_n(x)=\sum_{j=0}^n l_j(x)f(x_j)$
这个公式就是说根据输入的数据，我们得到了一个多项式。其实 $f(x_j)$ 是我们已经有的数据，也就是 $y_0,y_1,y_2,y_3$ 等等。
而复杂的就是 $l_j(x)$
$l_j(x)=\prod_{i=0,i\neq j}^n \frac{x-x_i}{x_j-x_i}=\frac{x-x_0}{x_j-x_0}\dots\frac{x-x_{j-1}}{x_j-x_{j-1}}\frac{x-x_{j+1}}{x_j-x_{j+1}}\dots\frac{x-x_{n}}{x_j-x_{n}}$
所以根据这个公式套着写代码就完了，我们先看 $l_j(x)$ 分母是一个常数，而分子很复杂，必须用数组表示，所以有以下类：
在这里插入图片描述
分式代码，分子是数组，分母是一个常数，Python代码如下：

# _*_ coding:utf-8 _*_

class LagrangeFraction:

    def __init__(self, numeractors:[], demominator):
        self.__numeractors = numeractors
        self.__demominator = demominator

    @property
    def numeractors(self):
        return self.__numeractors

    @property
    def demominator(self):
        return self.__demominator

每一项的代码：

# _*_ coding:utf-8 _*_

from com.youngthing.mathalgorithm.interpolation.lagrange_fraction import LagrangeFraction

# 拉格朗日多项式的一项
# 每一个分式分子是x - xi
# 每个分母是xj-xi,所有的分母可以乘起来，成为一个数字
# 所以是分子的list，分母是一个数字
class LagrangeItem:

    def __init__(self, y, lagrange_fraction:LagrangeFraction):
        # 有两个属性一个是y值 另一个是分式的连乘
        self.__y = y
        self.__lagrange_fraction = lagrange_fraction

    @property
    def y(self):
        return self.__y

    @property
    def lagrange_fraction(self):
        return self.__lagrange_fraction

多项式代码：

# _*_ coding:utf-8 _*_

from com.youngthing.mathalgorithm.interpolation.lagrange_item import LagrangeItem
from com.youngthing.mathalgorithm.interpolation.lagrange_fraction import LagrangeFraction


class LagrangePolynomial:

    def __init__(self, items: list[LagrangeItem]):
        self.__items: list[LagrangeItem] = items

    @staticmethod
    def calc(item: LagrangeFraction, x):
        # 连乘起来
        numerator = 1
        for numeractor in item.numeractors:
            numerator *= (x - numeractor)
        return numerator / item.demominator

    # 本身是个函数
    def __call__(self, *args, **kwargs):
        result = 0
        items = self.__items
        for item in items:
            result += LagrangePolynomial.calc(item.lagrange_fraction, args[0]) * item.y
        return result

创建方法的代码：

# _*_ coding:utf-8 _*_
from com.youngthing.mathalgorithm.interpolation.lagrange_item import LagrangeItem
from com.youngthing.mathalgorithm.interpolation.lagrange_fraction import LagrangeFraction
from com.youngthing.mathalgorithm.interpolation.lagrange_polynominal import LagrangePolynomial


def create_fraction(list, j):
    #先求分子的数组,就是x值
    numerators = []
    # 分母
    domominator = 1
    for i in range(0, len(list)):
        if i != j:
            xi = list[i][0]
            numerators.append(xi)
            xj = list[j][0]
            domominator = domominator*(xj - xi)
    return LagrangeFraction(numerators, domominator)

def interpolate(data:list):
    # 入参是二维数组
    items = [LagrangeItem(point[1], create_fraction(data, i)) for (i, point) in enumerate(data)]
    return LagrangePolynomial(items)

测试类：

if __name__ == '__main__':

    #测试数据
    data=[[0,1],[2,5],[4,17]]
    polynomial = lagrange.interpolate(data)
    print(polynomial(0),polynomial(2),polynomial(4))
    print(polynomial(1),polynomial(3),polynomial(5))

测试结果

1.0 5.0 17.0
2.0 10.0 26.0

然后再进行多项式化简，拉格朗日插值算法，获取的多项式，运算起来性能并不高。比如刚才的例子，得到的多项式是连乘的形式。写个简单的拉格朗日多项式的字符串方法：

def __str__(self):
    str = ''
    items = self.__items
    first = True
    for item in items:
        # 每项是一个(x-xn)的连乘*y/分母
        pi = ''
        for numeractor in item.lagrange_fraction.numeractors:
            pi += f'(x-{numeractor})'
        item_str = f'{pi}*{item.y}/({item.lagrange_fraction.demominator})'
        if (first) :
            first = False
        else:
            str += '+'
        str += item_str
    return str

上面的例子，得到的结果是这个样子：

(x-2)(x-4)*1/(8)+(x-0)(x-4)*5/(-4)+(x-0)(x-2)*17/(8)

所以，肯定要化简的啊，如果更高次的多项式，这个式子就更复杂了。我们第一步化简的是乘然后除这东西。对于这种多项式，还需要用到前面的知识对多项式进行化简，化简后成为这个样子：
$(x-2)\times(x-4)\times0.125+(x-0)\times(x-4)\times-1.25+(x-0)\times(x-2)\times2.125$
对于(x-2)*(x-4)*0.125这样的，可以使用多项式展开算法，可以变成下列的矩阵：
$\begin{pmatrix} 1 & -2\\ 1 & -4\\ 0 & 0.125 \end{pmatrix}$
但是为了简便计算，阔以把0.125直接乘到第一个多项式里所以,就是：
$\begin{pmatrix} 0.125& -0.25\\ 1 &-4 \end{pmatrix}$
因为前面写过多项式代数运算的代码，所以直接拿来用，代码变成了以下样子：

def __init__(self, items: list[LagrangeItem]):
    polynomials = []
    for item in items:
        y = item.y / item.lagrange_fraction.demominator
        # 进一步对多项式进行化简
        # 生产一个多项式矩阵
        matrix = [[1, -numeractor] for numeractor in item.lagrange_fraction.numeractors]
        matrix[0] = [y * i for i in matrix[0]]
        polynomials.append(polynomial_simplify(matrix))
    # 拉格朗日多项式的次数是项数-1，而多项式数组长度是次数+1，所以恰好等于项数
    merged_polynomial = [0] * len(items)
    for polynomial in polynomials:
        for i, coefficient in enumerate(polynomial):
            merged_polynomial[i] += coefficient
    self.__items: list[LagrangeItem] = items
    self.__polynomial = merged_polynomial

其中的ps是引入的多项式运算的库，如下

from com.youngthing.mathalgorithm.polynomial import polynomial_simplify

最终重载的call方法变成了这样简介（利用了霍纳法则）

def __call__(self, *args, **kwargs):
    result = 0
    for coefficient in self.polynomial:
        result = result * args[0] + coefficient
    return result