数值分析-多项式插值方法小结

小刀来啦

已于 2022-05-27 15:21:24 修改

阅读量4.4k

点赞数

分类专栏：数值分析笔记 python 数学原理文章标签：算法数据分析

于 2022-03-01 21:12:24 首次发布

本文链接：https://blog.csdn.net/weixin_40519529/article/details/123215263

版权

数值分析-多项式插值方法小结

前言
插值的应用与唯一性
Lagrange插值法和逐次线性插值
- 代码实现
逐次线性插值
- 代码实现
- 逼近复杂函数
Newton插值法
- 代码实践
差分与等距节点差分插值
- 代码实现
Hermite插值方法
写在最后

前言

最近在学数值分析，把一些算法实现一下

具体内容详见 《数值分析》李庆扬、王能超、易大义（华中科大出版社）

插值的应用与唯一性

插值的应用：

对于一些复杂函数，我们往往只能知道某些具体位置的函数值（实际中数据还有误差），当我们需要得到全体定义域上任一点的数值时，便可以采用插值的方法
可以用多项式插值函数来拟合复杂的数据点，简化函数表达。（多项式历来最被研究者喜欢，可积可导，性质很好！）

对于给定的n个点，所有的 $x_i \neq x_j$ ，即所有的采样点均不同，存在且唯一的次数不超过n次的多项式经过所有采样点，称之为插值多项式。（思路：组成的n阶线性方程组的系数矩阵行列式为范德蒙行列式，行列式不为零，所以存在唯一解。）

Lagrange插值法和逐次线性插值

Lagrange线性插值
- Lagrange插值的思想是不去求解各x幂次的系数，直接把函数值y作为系数确定下来，反过来确定线性空间的基的形式！
- 考虑两点的线性插值
  
  已知 $f(x_k)=y_k, f(x_{k+1})=y_{k+1}$ ，设两点的一阶插值多项式为 $L_1(x)$ ，有 $L_1(x_k)=y_k, L_1(x_{k+1})=y_{k+1}$ 。
  
  注意插值多项式一定经过对应的插值节点
  
  按照Lagrange插值的思想（由直线的点斜式和两点式可以简单推来），我们可以得到：
  
  $\Large L_1(x)=\frac{x-x_{k+1}}{x_k-x_{k+1}}y_k+\frac{x-x_{k}}{x_{k+1}-x_k}y_{k+1}$
  
  其实可以利用插值多项式的零点和抽样delta的性质来得到多点Lagrange插值多项式的形式：
  - 零点条件：对于 $l_k(x)$ ，在除第k个采样点的点上始终为0
  - 抽样性质：对于 $l_k(x)$ ，在第k个采样点的点上为1
- 给出 $l_k(x)$ 的一般形式：
  
  $\Large l_k(x)=\frac{(x-x_0)...(x-x_{k-1})(x-x_{k+1})...(x-x_n)}{(x_k-x_0)...(x_k-x_{k-1})(x_k-x_{k+1})...(x_k-x_n)} \,\,\,(k=0,1,...,n)$
- Lagrange插值多项式表达式：
  
  $\Large L_n(x)=\sum_{k=0}^ny_kl_k(x)=\sum_{i=0}^n y_i \prod_{j=0,j \neq i}^n \frac{x-x_j}{x_i-x_j}$
  
  其中，我们可以看到Lagrange插值多项式的基函数为：
  
  $base_n(x)=l_k(x)=\prod_{j=0,j \neq k}^n \frac{x-x_j}{x_k-x_j}$
- 插值余项
  定义截断误差（插值余项）为 $R_n(x)=f(x)-L_n(x)$
  根据罗尔定理，可以推得：
  
  $\large R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^n (x-x_i)$

代码实现

import numpy as np
from functools import reduce

# Lagrange 线性插值函数
class Lagrange_Polynomial_Interpolation(object):
    """Lagrange_Polynomial_Interpolation 

    base function: Ln(x)=\sum_{k=0}^{n} yk l_k(x)
    Author: Junno
    Date: 2022-03-01
    """

    def __init__(self, x, y):
        """__init__ 

        Args:
            x ([np.array]): [x of points]
            y ([np.array]): [y of points]
        """
        self.N = x.shape[0]
        self.x = np.copy(x)
        self.y = np.copy(y)

    def calculate_lk(self, x, k,):
        """calculate_lk(x)

        base function: lk(x)=∏——{j=0,j!=k}^{n} (x-xj)/(xk-xj)
        Args:
            x ([float]): [x]
            k ([int]): [k]

        Returns:
            [float]: [lk(x)]
        """
        ta = reduce(lambda x, y: x*y, [x-self.x[i]
                    for i in range(self.N) if i != k])
        tb = reduce(lambda x, y: x*y, [self.x[k]-self.x[i]
                    for i in range(self.N) if i != k])
        return ta/tb

    def predict(self, x):
        """predict 

        Args:
            x ([float]): [x]

        Returns:
            [float]: [predict y]
        """
        return np.sum([self.y[k]*self.calculate_lk(x, k) for k in range(self.N)])

## 测试教材例子
# f(x)=sin(x)
x=np.array([0.32,0.34,0.36])
y=np.array([0.314567,0.333487, 0.352274])

LPI=Lagrange_Polynomial_Interpolation(x,y)
x0=0.3367
LPI.predict(x0)
>>> 0.3303743620375
np.sin(x0)-LPI.predict(x0)
>>> -1.7048187195278786e-07

逐次线性插值

Lagrange插值多项式计算函数近似值时，如果加入了新的采样点，则所有权重表达式都需要重新计算，不是很高效，逐次线性插值法可以在已有计算的数据上得到高阶的权重表达，更加方便和高效。

一般情况下，两个k次的插值多项式可以通过线性插值得到一个k+1次的插值多项式 （具体推导详见教材21页）。

记 $I_{i_1,i_2,...,i_n}(x)$ 表示关于点 $x_{i_1},x_{i_2},...,x_{i_n}$ 的 $(n - 1)$ 次插值多项式 （n个数据点最多组成n-1次），特别地， $I_{i_k}(x)$ 为零次插值多项式，即函数值 $f(x_{i_k})$ 。

Neville算法：

$\Large I_{0,1,...,k+1}(x)=I_{0,1,...,k}(x)+\frac{I_{1,...,k+1}(x)-I_{0,1,...,k}}{(x_{k+1}-x_0)}(x-x_0)$

可以理解为要计算0到k+1点的k+1次插值多项式时，可以由0到k点的k次插值多项式和1到k+1点的k次插值多项式通过线性插值得到，这样便组成了一个可以递归求解的公式。

代码实现

# 逐次线性插值
class mini_I(object):
    """
    class for calculating the interpolation of order k+1 from two component of order k
    base funtion: ans=y0+(y1-y0)/(x1-x0)*(x-x0)

    Author: Junno
    Date: 2022-03-01
    """

    def __init__(self, x0=0, x1=0, y0=0, y1=0, first=False):
        """__init__ 
        Args:
            x0 (float, optional): [component of the base function]. Defaults to 0.
            x1 (float, optional): [component of the base function]. Defaults to 0.
            y0 (float, optional): [component of the base function]. Defaults to 0.
            y1 (float, optional): [component of the base function]. Defaults to 0.
            first (bool, optional): [whether the first element of each order list]. Defaults to False.
        """
        self.x0 = x0
        self.x1 = x1
        self.y0 = y0
        self.y1 = y1
        self.first = first
        # record answer for fast calculation
        self.last_answer = {
   }

    def __call__(self, x=0.):
        """__call__ 

        Args:
            x ([float]): [interpolation node]

        Returns:
            [answer]: [recursive call or the recorded values]

        """
        if x not in self.last_answer.keys():
            # here x0,x1,y0,y1 can be a mini_I class to call for input x
            answer = self.y0(x)+(self.y1(x)-self.y0(x))/(self.x1 -
                                                         self.x0)*(x-self.x0) if not self.first else self.y0
            # record
            self.last_answer[x] = answer
            return answer
        else:
            if x in self.last_answer.keys():
                return self.last_answer[x]


class Successive_Linear_Interpolation_2D(object):
    """Successive_Linear_Interpolation_2D 逐次线性插值

    Author: Junno
    Date: 2022-03-01
    """

    def __init__(self, x, y, base_order=0):
        """__init__ 

        Args:
            x ([np.array]): [x of points]
            y ([np.array]): [y of points]
            base_order ([int]): [interpolation order]
        """
        self.N = x.shape[0]
        assert base_order < self.N, 'The base_order should be less than the number of points'
        self.x = np.copy(x)
        self.y = np.copy(y)
        self.base_order = 0
        # flag of calculated
        self.calculated = False

        # ref:《数值分析》 李庆扬、王能超、易大义（华中科大出版社）
        self.F = lambda i, j: mini_I(
            self.x[i-j], self.x[i], self.I_set[i-1][j-1], self.I_set[i][j-1], first=False)

    def calculate_I(self, x, eps=1e-6, order=1):
        """calculate_I calculate y for a given x and order with threshold eps

        Args:
            x ([float]): [x]
            eps ([float], optional): [threshold to stop process]. Defaults to 1e-6.
            order ([int], optional): [interpolation order]. Defaults to 1.
        """
        assert order < self.N, 'The order should be less than the number of points'

        if not self.calculated:

            self.base_order = self.N if order is None else order

            self.I_set = [[]]

            self.I_set[0].append(mini_I(y0=self.y[0], first=True))

            last = 0

            for i in range(1, self.base_order+1):
                self.I_set += [[]]
                self.I_set[i].append(mini_I(y0=self.y[i], first=True))
                for k in range(1, i+1):
                    self.I_set[i].append(self.F(i, k))

                # early stop when the difference between of two results from k order and k+1 order is less than eps
                if i > 1:
                    new = self.I_set[i][-1](x)
                    # print(i, abs(last-new) <= eps)
                    if abs(last-new) <= eps:
                        if i < (self.base_order+1)-1:
                            print(
                                "Meet with absotion condition of eps={} and the current order of interpolation is {}".format(eps, i))
                            self.base_order = i
#                             self.I_set.pop(-1)

                        self.calculated = True
                        return new

                last = self.I_set[i][-1](x)

                # print(i,)
                # self.print(x)

            self.calculated = True

            return self.I_set[self.base_order][self.base_order](x)

        else:
        	last = self.I_set[0][0](x)
            for i in range(1, order+1):

                if i > self.base_order:
                    self.I_set += [[]]
                    self.I_set[i].append(mini_I(y0=self.y[i], first=True))
                    for k in range(1, i+1):
                        self.I_set[i].append(self.F(i, k))

                new = self.I_set[i][-1](x)
                if abs(last-new) <= eps:
                    if check:
                        print(
                            "Meet with absotion condition of eps={} and the current order of interpolation is {}".format(eps, i))
                    return new
                last = new

            return self.I_set[order][order](x)

    def add_point(self, x, y):
        """add_point 

        Args:
            x ([np.array]): [x of points]
            y ([np.array]): [y of points]
        """
        M = x.shape[0]
        self.x = np.concatenate((self.x, x), axis=0)
        self.y = np.concatenate((self.y, y), axis=0)
        self.N = self.N+M

    def print(self, x):
        """print show I_set for a given x

        Args:
            x ([float]): [x]
        """
        for i in range(len(self.I_set)):
            print('I_{}:'.format(i), end='\t')
            for j in range(len(self.I_set[i])):
                print(self.I_set[i][j](x), end='\t')
            print('\n')

## 测试样例
import numpy as np
# f(x)=sin(x)
x=np.array([0.32,0.34,0.36])
y=np.array([0.314567,0.333487, 0.352274])

SLI=Successive_Linear_Interpolation_2D(x,y)

SLI.calculate_I(0.3345, order=2)
>>> 0.32829725843749996

SLI.calculate_I(0.3345, order=2)-np.sin(0.3345)
>>> 3.3479488331655816e-07

# print I_set structure
SLI.print(0.3345)
>>>
I_0:	0.314567	

I_1:	0.333487	0.32828399999999996	

I_2:	0.352274	0.32832057499999995	0.32829725843749996	


# add points
SLI.add_point(np.array([0.38+i*0.02 for i in range(3)]),np.array([np.sin(0.38+i*0.02) for i in range(3)]))

SLI.x
>>> array([0.32, 0.34, 0.36, 0.38, 0.4 , 0.42])

# calculate x=0.3678 of order 5-th
SLI.calculate_I(0.3678, order=5)
>>> 
Meet with stopping condition of eps=1e-06 and the current order of interpolation is 4
0.3595632718504201

SLI.print(0.3678)
>>>
I_0:	0.314567	

I_1:	0.333487	0.35978579999999993	

I_2:	0.352274	0.35960093000000004	0.35956488035000006	

I_3:	0.3709204694129827	0.35954612307106326	0.35956283918438897	0.35956325422139657	

I_4:	0.3894183423086505	0.35963676694662533	0.3595637986267979	0.3595632837260384	0.3595632718504201

逼近复杂函数

# 使用插值多项式逼近复杂函数
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
x=np.arange(0,1,0.01) # x.shape=200
F=lambda x:np.cos(2*np.pi*x)+np.sin(x)+np.exp(x**2) # original_function

# sample x
xs=np.sort(np.random.choice(x, size=(10), replace=False))

# with/without measurement error
ys=F(xs)
ys_noise=F(xs)+np.random.randn(xs.shape[0],)*0.05

SLI=Successive_Linear_Interpolation(xs,ys)
SLI_n=Successive_Linear_Interpolation(xs,ys_noise)

y_n=np.zeros_like(x)
y=np.zeros_like(x)

for k in range(x.shape[0]):
    y[k]=SLI.calculate_I(x[k],order=xs.shape[0]-1,check=False)
    y_n[k]=SLI_n.calculate_I(x[k],order=xs.shape[0]-1,check=False)

# show results
plt.figure(figsize=(8,5))
S=plt.scatter(xs,ys,s=5, marker='*',c='black')
Predict_without_noise=plt.plot(x,y,c='red')
Predict_with_noise=plt.plot(x,y_n,c='yellow')
# Raw=plt.plot(x,F(x), c='g')
plt.xlabel('x')
plt.ylabel('y')
Samples=mpatches.Patch(color='black', label='Samples Points')
Predict_without_noise=mpatches.Patch(color='red', label='Predict_without_noise')
Predict_with_noise=mpatches.Patch(color='yellow', label='Predict_with_noise')
plt.title('Successive_Linear_Interpolation with/withour noise')
plt.legend(handles=[Samples,Predict_without_noise,Predict_with_noise],loc='best')
plt.show()

在这里插入图片描述

可以看到线性插值对噪声还是非常敏感的，可以保证的是插值多项式一定过采样点，但是由于噪声两边会剧烈变化。可见，线性插值方法只有在插值节点数据足够准确时才可以比较好地逼近原函数。
当加入的噪声太大时，有可能出现以下的情况：

在这里插入图片描述

Newton插值法

本质上由于只有n+1个点的情况下，线性插值多项式存在且唯一存在。 不同插值方法，如Lagrange、Neville、Newton等都是从不同的基函数入手来对插值多项式进行化简总结，本质上得到的插值多项式以及插值余项在数学上是等价的。在Lagrange的基础上，Neville和Newton插值法可以利用已计算的阶数数据来求解下一阶的内容，更加高效便捷。
一般Newton插值法由差商概念引入，这里从最本质的Lagrange插值法推出Newton插值法

首先，Lagrange插值法的n+1点插值多项式表达式如下：

$\Large L_n(x)=\sum_{k=0}^ny_kl_k(x)=\sum_{i=0}^n y_i \prod_{j=0,j \neq i}^n \frac{x-x_j}{x_i-x_j}$

最低0.47元/天解锁文章