拉格朗日插值与牛顿插值

Posted 2020-11-25 philolif

拉格朗日插值与牛顿插值

这里使用python实现了拉格朗日插值和牛顿插值

import numpy as np


# 拉格朗日插值
class Lagrange:
    def __init__(self, x, y):
        """
        这里x,y是长度相同的一维numpy数组
        """
        # 插值基函数的系数
        self.coefficient = []
        # 基函数
        self.base = []
        for i in range(x.size):
            self.base.append(np.delete(x, i))
            self.coefficient.append(
                y[i] / np.prod(x[i] - self.base[i])
            )
        self.base = np.array(self.base)
        self.coefficient = np.array(self.coefficient)

    def calculate(self, x):
        """
        这里的x是一个标量
        """
        return np.sum(
            np.prod(x - self.base, axis=1) * self.coefficient
        )


# 牛顿插值
class Newton:
    def __init__(self, x, y):
        """
        这里x,y是长度相同的一维numpy数组
        """
        # 插值节点的个数
        self.n = x.size
        # 差商表
        self.table = np.zeros((self.n, self.n + 1))
        self.table[:, 0] = x
        self.table[:, 1] = y
        # 按行构造差商表
        for i in range(1, self.n):
            for j in range(2, i + 2):
                self.table[i, j] = (self.table[i, j - 1] - self.table[i - 1, j - 1]) / (
                        self.table[i, 0] - self.table[i - j + 1, 0])

    def calculate(self, x):
        """
        这里x是一个标量
        """
        value = 0
        for i in range(1, self.n):
            value += self.table[i, i + 1] * np.prod(x - self.table[:i, 0])
        # 加上f(x0)
        value += self.table[0, 1]
        return value

    def append(self, x, y):
        """
        牛顿插值可以插入更多的点以提高精度
        这里x,y是长度相同的一维numpy数组
        """
        n = self.n + x.size
        table = np.zeros((n, n + 1))
        # 对原有的差商表进行扩大
        table[:self.n, :self.n + 1] = self.table
        table[self.n:, 0] = x
        table[self.n:, 1] = y
        # 按行构造差商表
        for i in range(self.n, n):
            for j in range(2, i + 2):
                table[i, j] = (table[i, j - 1] - table[i - 1, j - 1]) / (table[i, 0] - table[i - j + 1, 0])
        self.n = n
        self.table = table


if __name__ == '__main__':
    # example
    # 拉格朗日插值
    lx = Lagrange(
        x=np.array([1, 3, 4, 7, 9], dtype=np.float64),
        y=np.array([0, 2, 15, 12, 2], dtype=np.float64)
    )
    print("L(2.5)=", lx.calculate(2.5))
    # 3个点进行牛顿插值
    nx = Newton(
        x=np.array([1, 3, 4], dtype=np.float64),
        y=np.array([0, 2, 15], dtype=np.float64)
    )
    print("N(2.5)=", nx.calculate(2.5))
    # 牛顿插值新增2个点
    nx.append(
        x=np.array([7, 9], dtype=np.float64),
        y=np.array([12, 2], dtype=np.float64)
    )
    print("N(2.5)=", nx.calculate(2.5))

测试结果

L(2.5)= -3.98203125
N(2.5)= -1.5
# 新增点后的牛顿插值结果和拉格朗日插值结果相同
N(2.5)= -3.98203125