数值计算系列之牛顿插值原理及实现

原创

devecor

修改于 2020-05-25 10:48:51

2.5K0

修改于 2020-05-25 10:48:51

文章被收录于专栏：数学与计算机

前言

本文介绍了插值概念和一般的计算方法，介绍了用于简化插值函数计算的牛顿插值方法。最后给出牛顿向前插值算法的python实现。

若您在阅读的过程中发现任何问题, 邮箱联系我: devecor@163.com

所谓插值？

假定存在一个未知函数 y=f(x) ,已知若干个点处的函数值：

f(x)满足上述数据，但是满足上述数据的函数有无数个，没有办法从这一组数据中准确的得出f(x)

所谓插值，就是找到一个满足上述数据的足够好的函数，来代表或者近似f(x)，笔者将这样的一个插值函数记为\phi(x)

什么是足够好的函数？

一个足够好的近似函数\phi(x)，应当具备以下特点：

形式简单
容易计算
足够好的性态，即具有足够高阶的导数（连续优于间断，可导优于连续）

具有这样的特点的函数，通常是一类函数，比如多项式。

计算插值函数

确定一个函数类作为未知函数f(x)的近似函数。这个函数类可以是多项式，可以是三角函数，也可以是任何你认为具有上述特点的某一类函数。
从这个函数类中找到一个满足\phi(x_i)=y_i(i=1,2,...n)条件的函数\phi(x)作为这组数据的插值函数。前述条件称为插值条件，这些数据称为插值节点

多项式插值

一旦我们选定了一个函数类，就能知道所求的插值函数具有怎样的形式。比如选定多项式作为近似函数，使用下面的数据：

该组数据一共有n+1个插值节点，也就是说具有n+1个插值条件，要使这些插值条件同时得到满足，只能使用n次多项式作为插值函数，显然这样的插值函数具有以下形式：

y=a_0+a_1x+a_2x^2+\cdots +a_nx^n

根据插值条件，接下列方程即可

容易证明只要插值节点互不相同，上述方程必然有唯一解。

那么，问题来了，你愿不愿意求解上述的方程组呢？假如有2000条数据，即有2000个方程的高次方程组，你解过这么大的吗？或者说这种方法是有效的吗？

简化计算：牛顿插值

假定只有一个插值节点，那么数据看起来像这样

只有一个数据的插值函数记为\phi_0(x),要满足插值条件，插值函数显然是

\phi_0(x)=y_0

增加一个插值节点，数据像这样：

插值函数记为\phi_1(x)可以写成下面的形式：

\phi_2(x)=y_0+a_1(x-x_0)

上式满足第一个插值节点，a_1未知，将第二个插值节点代入，求得

a_1=\frac{y_1-y_0}{x_1-x_0}=\frac{y_1-\phi_0(x_1)}{x_1-x_0}

再增加一个插值节点，数据：

插值函数\phi_3(x)=y_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)

上式满足前两个插值节点，a_2未知，代入第三个插值节点，求得a_2=\frac{y_2-\left[y_0+a_1(x_2-x_0)\right]}{(x_2-x_0)(x_2-x_1)}=\frac{y_2-\phi_1(x_2)}{(x_2-x_0)(x_2-x_1)}

一般的，

a_i=\frac{y_i-\phi_{i-1}(x_i)}{(x-x_0)(x-x_1)\cdots(x-x_{i-1})}

那么，

\phi_n(x) = y_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+\cdots+a_n(x-x_0)(x-x_1)\cdots(x-x_{n-1})

非等距节点牛顿插值的python实现

若插值节点是等距的，则可以使用向前差分极大的简化运算。不失一般性，本文基于递归实现非等距节点的牛顿插值算法。

功能设计

Newton插值算法极大的简化了计算，相对于Lagrange插值，最大的优势在于其具有计算的继承性。所谓继承性是指当增加一条数据，不必与原有数据一起重新计算插值函数，而是在原有的插值函数上增加一项，得出新的插值函数。基于此我们的算法应有下面的几个功能：

数据的输入
数据的更新
计算插值函数，保存计算结果

所以，新建pyiplt.py文件，写下如下代码

class Interpolate(object):
    """Newton 插值函数类"""
    def __init__(self, data):
        pass
    
    def update(self, data):
        pass
    
    def interpolate(self):
        """计算插值函数"""
        pass

数据的输入和存在形式

首先需要定义数据输入为nx2的numpy数组，为使用方便对__init__()函数做如下改动, 并定义数据的getter:

class Interpolate(object):

    def __init__(self, data=None, x=None, y=None):
        pass

    @property
    def data(self):
        return self.__data

新建 test_pyiplt.py文件

测试输入数据存在形式, 写下单元测试：

from unittest import TestCase
import numpy as np
from pyiplt import Interpolate


class TestInterpolate(TestCase):
    x = [i for i in range(10)]
    y = x.copy()
    res = np.array([[i, i] for i in range(10)]).reshape((-1, 2))

    def test_data(self):

        interpolator = Interpolate()
        self.assertIsNone(interpolator.data)

        interpolator = Interpolate(x=self.x, y=self.y)
        self.assertEqual(self.res.tolist(), interpolator.data.tolist())

刚好能通过测试的代码：

import numpy as np

class Interpolate(object):

    def __init__(self, data=None, x=None, y=None):

        if data is not None:
            self.__data = data
        elif x is not None and y is not None:
            assert len(x) == len(y), "lengths isn't equal"
            self.__data = np.array([[i, j] for i, j in zip(x, y)])\
                .reshape((-1, 2))
        else:
            self.__data = None
    
    @property
    def data(self):
        return self.__data

计算插值函数

先写单元测试：

from unittest import TestCase
from pyiplt import Interpolate


class TestInterpolate(TestCase):

    x = [i for i in range(10)]

    y2 = [f(i) for i in x]
    y3 = [f2(i) for i in x]

    res2 = [1, 2, ] + [0 for _ in range(8)]    # 升幂排列
    res3 = [2, 8, 3] + [0 for _ in range(7)]   # 升幂排列

    def test_interpolate(self):

        interpolator = Interpolate()
        self.assertRaises(Exception, interpolator.interpolate)

        interpolator = Interpolate(x=self.x, y=self.y2)
        self.assertEqual(self.res2, interpolator.interpolate())

        interpolator2 = Interpolate(x=self.x, y=self.y3)
        self.assertEqual(self.res3, interpolator2.interpolate())

在给出能够通过测试的代码之前，交代以下实现的思路：

根据Newton插值公式：

\phi_n(x) = y_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+\cdots+a_n(x-x_0)(x-x_1)\cdots(x-x_{n-1}) a_i=\frac{y_i-\phi_{i-1}(x_i)}{(x-x_0)(x-x_1)\cdots(x-x_{i-1})}

我们知道，只需计算出a_0, a_1, ..., a_n就可以得到插值函数，但是a_n与前n个数据（n+1条数据）的插值函数有关，也就是与前面n个系数有关（你品，你细品），直到a_0，我们知道a_0 = y_0，只跟第一个插值节点相关。解决这种问题最常见的思维方式是自顶向下的模式，计算n+1个数据的插值函数时，先去计算前n个数据的插值函数，先计算前n-1个数据的插值函数，直到最后只有一个插值节点，写出这个节点的插值函数，然后依次回溯。这样的方式叫递归，通过对自身的调用化繁为简，最后逐步回溯。

import numpy as np


def sub_mult(data):
    """
    计算系数a_n分母处的连乘
    :param data: 插值节点序列, [x0, x1, x2, x3, ..., x_n]
    :return: (x_n - x_0) * (x_n - x_1) * ... * (x_n - x_{n-1})
    """
    res = 1
    loop_times = len(data) - 1
    i = 0
    while i < loop_times:
        res *= (data[-1] - data[i])
        i += 1
    return res


class Interpolate(object):
    """Newton插值函数类
    非等距节点的Newton插值的递归实现
    :param data: array_like: [[x1, y1], [x2, y2], ...]
    :param x: 插值节点序列
    :param y: 对应函数值序列
    """

    def __init__(self, data=None, x=None, y=None):
        if data is not None:
            self.__data = data
        elif x is not None and y is not None:
            assert len(x) == len(y), "lengths isn't equal"
            self.__data = np.array([[i, j] for i, j in zip(x, y)])\
                .reshape((-1, 2))
        else:
            self.__data = None

        # 插值结果的系数列表
        self.__coefficients = []

    def interpolate(self):
        """
        计算插值函数, 返回系数列表
        :return: list 多项式插值函数系数列表，按升幂排列
        """
        if self.__data is None or len(self.__data) == 0:
            raise Exception("NoneDataError: please update first, use 'self.update'")
        self.__factor(self.__data)
        return self.__coefficients

    def __factor(self, data):
        """
        自顶向下递归的计算插值函数的系数
        :param data:
        :return:
        """
        x = [i[0] for i in data]
        if data.shape[0] == 1:
            self.__coefficients.append(data[0][1])
            return
        self.__factor(data[0:-1])
        self.__coefficients.append(
            (data[-1][1] - self.__phi(x)) / sub_mult(x)
        )

    def __phi(self, data):
        """
        计算f_{n-1}(x_n)的值
        :param data: [x_0, x_1, ..., x_n]
        :return:
        """
        assert len(data) >= 2, "assertFail: length < 2!!!"
        res = self.__coefficients[0]
        if len(data) > 2:
            loop_times = len(data[0:-2])
            i = 0
            while i < loop_times:
                temp = 1
                j = 0
                while j <= i:
                    temp *= (data[-1] - data[j])
                    j += 1

                res += self.__coefficients[i+1] * temp
                i += 1
        return res

下面是完整版的代码，比前文更完善，请读者细品

单元测试：

import numpy as np
from unittest import TestCase
from pyiplt import Interpolate
from pyiplt import sub_mult


def f(x):
    return 2 * x + 1


def f2(x):
    return 2 + 5 * x + 3 * x ** 2


class TestInterpolate(TestCase):

    x = [i for i in range(10)]
    y = x.copy()

    y2 = [f(i) for i in x]
    y3 = [f2(i) for i in x]

    res = np.array([[i, i] for i in range(10)]).reshape((-1, 2))
    res2 = [1, 2, ] + [0 for _ in range(8)]    # 升幂排列
    res3 = [2, 8, 3] + [0 for _ in range(7)]   # 升幂排列

    def test_data(self):

        interpolator = Interpolate()
        self.assertIsNone(interpolator.data)

        interpolator = Interpolate(x=self.x, y=self.y)
        self.assertEqual(self.res.tolist(), interpolator.data.tolist())

    def test_interpolate(self):

        interpolator = Interpolate()
        self.assertRaises(Exception, interpolator.interpolate)

        interpolator = Interpolate(x=self.x, y=self.y2)
        self.assertEqual(self.res2, interpolator.interpolate())

        interpolator2 = Interpolate(x=self.x, y=self.y3)
        self.assertEqual(self.res3, interpolator2.interpolate())

    def test_iplt_f(self):
        interpolator = Interpolate(x=self.x, y=self.y2)
        interpolator.interpolate()
        res = [interpolator.iplt_f(i) for i in self.x]
        self.assertEqual(res, self.y2)


class Test(TestCase):
    def test_sub_mult(self):
        data1 = [1, 2, 3]
        res1 = 2
        data2 = [3, 6, 4]
        res2 = -2
        data3 = [8, 3, 5, 9]
        res3 = 24
        self.assertEqual(res1, sub_mult(data1))
        self.assertEqual(res2, sub_mult(data2))
        self.assertEqual(res3, sub_mult(data3))

可以通过测试的代码：

import numpy as np


def sub_mult(data):
    """
    计算系数a_n分母处的连乘
    :param data: 插值节点序列, [x0, x1, x2, x3, ..., x_n]
    :return: (x_n - x_0) * (x_n - x_1) * ... * (x_n - x_{n-1})
    """
    res = 1
    loop_times = len(data) - 1
    i = 0
    while i < loop_times:
        res *= (data[-1] - data[i])
        i += 1
    return res


class Interpolate(object):
    """Newton插值函数类
    非等距节点的Newton插值的递归实现
    :param data: array_like: [[x1, y1], [x2, y2], ...]
    :param x: 插值节点序列
    :param y: 对应函数值序列
    """

    def __init__(self, data=None, x=None, y=None):
        if data is not None:
            self.__data = data
        elif x is not None and y is not None:
            assert len(x) == len(y), "lengths isn't equal"
            self.__data = np.array([[i, j] for i, j in zip(x, y)])\
                .reshape((-1, 2))
        else:
            self.__data = None

        # 插值结果的系数列表
        self.__coefficients = []

    @property
    def data(self):
        return self.__data

    def update(self, data=None, x=None, y=None):
        """
        在现有数据上添加数据，返回更新后的插值结果，充分发挥Newton插值的优势
        :param data: array_like: [[x1, y1], [x2, y2], ...]
        :param x: 插值节点序列
        :param y: 对应函数值序列
        :return: list 多项式插值函数系数列表，按升幂排列
        """

        # 拼接数据
        if data is not None:
            self.__data = np.concatenate((self.__data, data))
        elif x is not None and y is not None:
            assert len(x) == len(y), "lengths isn't equal"
            data = np.array([[i, j] for i, j in zip(x, y)])\
                .reshape((-1, 2))
            self.__data = np.concatenate((self.__data, data))

        # 更新插值函数
        i = data.shape[0]
        while i > 0:
            if i-1 == 0:
                x = [i[0] for i in self.__data]
            else:
                x = [i[0] for i in self.__data[0:-(i-1)]]

            self.__coefficients.append(
                (self.__data[-(i-2)][1] - self.__phi(x)) / sub_mult(x)
            )
        return self.__coefficients

    def interpolate(self):
        """
        计算插值函数, 返回系数列表
        :return: list 多项式插值函数系数列表，按升幂排列
        """
        if self.__data is None or len(self.__data) == 0:
            raise Exception("NoneDataError: please update first, use 'self.update'")
        self.__factor(self.__data)
        return self.__coefficients

    def __factor(self, data):
        """
        自底向上递归的计算插值函数的系数
        :param data:
        :return:
        """
        x = [i[0] for i in data]
        if data.shape[0] == 1:
            self.__coefficients.append(data[0][1])
            return
        self.__factor(data[0:-1])
        self.__coefficients.append(
            (data[-1][1] - self.__phi(x)) / sub_mult(x)
        )

    def __phi(self, data):
        """
        计算f_{n-1}(x_n)的值
        :param data: [x_0, x_1, ..., x_n]
        :return:
        """
        assert len(data) >= 2, "assertFail: length < 2!!!"
        res = self.__coefficients[0]
        if len(data) > 2:
            loop_times = len(data[0:-2])
            i = 0
            while i < loop_times:
                temp = 1
                j = 0
                while j <= i:
                    temp *= (data[-1] - data[j])
                    j += 1

                res += self.__coefficients[i+1] * temp
                i += 1
        return res

    def iplt_f(self, x):
        """计算好的插值函数"""
        if len(self.__coefficients) > 0:
            res = 0
            for k, v in enumerate(self.__coefficients):
                i = 0
                temp = 1
                while i < k:
                    temp *= (x - self.__data[i][0])
                    i += 1
                res += v * temp
            return res
        else:
            raise Exception("UserError: please call self.interpolate first")

    def res2str(self):
        res = ""
        for i, e in enumerate(self.__coefficients):
            if i == 0:
                res += " " + str(e) + "\n"
            elif e > 0:
                res += ("+ " + str(e) + "x^" + str(i) + "\n")
            elif e < 0:
                res += ("- " + str(-e) + "x^" + str(i) + "\n")
        return res


if __name__ == "__main__":

    def f(x):
        return np.power(x, 2.6)

    x = np.array([i/10 for i in range(20, 31)])
    y = np.array([f(i) for i in x])

    print("输入数据：")
    print("x: {}".format(x))
    print("y: {}".format(y))

    interpolator = Interpolate(x=x, y=y)
    coefficients = interpolator.interpolate()

    print("插值函数:\n{}".format(interpolator.res2str()))

    print("计算结果：")
    print("2.34^2.6 = {}".format(interpolator.iplt_f(2.34)))
    print("2.98^2.6 = {}".format(interpolator.iplt_f(2.98)))

单元测试结果：