【数值分析】第2章-函数插值方法及其Python实现

1、前言

这个学期，上了一门叫做《应用数值分析(华南理工大学出版社)》(如下图所示)的数学课，为了巩固一下所学知识，同时也为了给其他同学一个参考，特撰写《数值分析》系列博客，主要包括：

【数值分析】第2章-函数插值方法及其Python实现
【数值分析】第3章-曲线拟合/函数逼近及其Python实现
【数值分析】第4章-数值微分/数值积分及其Python实现
【数值分析】第5章-线性代数方程组数值解法——直接法及其Python实现
【数值分析】第6章-线性代数方程组数值解法——迭代法及其Python实现
【数值分析】第7章-非线性代数方程组数值解法及其Python实现
【数值分析】第9章-常微分方程数值解法及其Python实现

本篇主要讲解函数插值方法然后使用Python实现所讲解的算法。

2、编程环境介绍及环境搭建

Python由于其拥有大量的开源包(例如numpy、scipy、matplotlib)，最近几年，颇受各领域的科研工作者追捧。本系列文章都是使用python编程实现的。
受篇幅影响，本文不详细介绍如何安装搭建Python编程环境，对这部分毫无经验的同学可以参考我之前撰写的博客，使用Windows的同学：【Python学习】Windows10开始你的Anaconda安装与Python环境管理、使用Linux的同学：【Python学习】纯终端命令开始你的Anaconda安装与Python环境管理，使用Mac的同学则更方便，直接在anaconda官网下载相应的pkg，然后安装即可。

3、插值问题的提出

科学背景

在科学研究或者工程问题当中，我们通常会遇到下列情况：

• 已知函数太复杂，难以输入到计算机，或者利用计算机运行一次，需要进行大量计算；
• 已知的函数只是一张离散的表格，也就是一些$(x_i,y_i)$点(由实验或观测得到)；
• 已知的数据往往会存在不完整的情况，希望通过某种方法补全缺失的数据。

插值和拟合的区别

这两种情况下，我们通常希望寻求简单的函数来逼近原函数。那么，有的同学可能就会问：这和拟合有什么区别？如果不知道什么是拟合，可以跳过这一小节。

拟合长这个样子：

插值长这个样子：

更直白地说，拟合通常不要求拟合函数完全穿过所给的数据点，但是插值，则通常要求插值函数经过所给的数据点。

4、插值问题的数学知识

插值问题

定义：设$f$是定义在区间$[a,b]$上的实值函数，并已知在$[a,b]$上$n+1$个互异节点$x_i$及相应的函数值$f_i=f(x_i)(i=0,1,\cdots ,n)$，要求建造一个简单的、便于计算的满足下列条件的函数$p$：
$$p(x_i)=f_i(i=0,1,\cdots ,n)$$

从而以$p$作为$f$的近似函数，其中
$f$：被插函数；
$p$：插值函数；
$x_i$：插值节点；
$[a,b]$：插值区间。

多项式插值的存在唯一性定理

设已知$[a,b]$上的函数$f$在$n+1$个互异节点$x_i\in [a,b]$上的值$f_i=f(x_i)(i=0,1,\cdots ,n)$，则存在唯一的次数$\le n$的多项式$p(x)\in P_n$满足：
$$p_n(x_i)=f_i(i=0,1,\cdots ,n)$$

5、常见插值公式

5.1 Lagrange插值公式

5.1.1 插值基函数与插值函数

根据多项式插值的存在唯一性定理，已知$[a,b]$上的函数$f$在$n+1$个互异节点$x_i\in [a,b]$，则必存在唯一的次数$\le n$的多项式$L_n(x)$满足：
$$L_n(x_i)=f_i(i=0,1,\cdots ,n)$$
定义插值基函数$l_i(x)$：
$$l_i(x)=\frac{(x-x_0)\cdots (x-x_{i-1})\cdots (x-x_n)}{(x_i-x_0)\cdots (x_i-x_{i-1})\cdots (x_i-x_n)}=\prod _{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}(i=0,1,\cdots ,n)$$
则插值多项式可以写为：
$$L_n(x)=\sum _{i=0}^n{l}_i(x)\cdot f_i$$
相对应的插值余项(用于误差估计)：
$$R_n(X)=\frac{1}{(n+1)!}{f}^{(n+1)}(\xi ){\omega }_{n+1}(x)$$
其中，${\omega }_{n+1}(x)={\prod }_{j=0}^{n}x-x_j$

5.1.2 Lagrange插值的Python实现

1、定义插值基函数

class DimException(Exception):
    def __init__(self, X,Y):
        err = 'X和Y的维度不相等,请检查数据输入是否正确!'
        Exception.__init__(self, err)
        self.X = X
        self.Y = Y


def lagrange_inter_base(X,Y,x):
    '''
    n次Lagrange插值的基函数数值计算
    X:已知的n+1个点的x_i
    Y:已知的n+1个点的y_i
    X和Y的维度必须一致，否则报错
    x:待通过插值函数L(x)计算具体函数值的x
    Return:
    l:插值基函数计算x后返回的结果，它是一个列表，维度和X、Y一致
    '''
    if len(X)!=len(Y):
        raise DimException(X,Y)
    
    n = len(Y)
    l = []
    for i in range(n):
        li = 1 #基函数初始值设置为1，方便进行循环求乘
        # 计算第i个插值基函数的值
        for j in range(n):
            if j!=i:
                li=li*(x-X[j])/(X[i]-X[j])
        l.append(li)
    return l

2、定义插值函数

def lagrange_inter(X,Y,x):
    '''
    Lagrange插值
    n次Lagrange插值的数值计算
    X:已知的n+1个点的x_i
    Y:已知的n+1个点的y_i
    x:一维变量
    '''
    if len(X)!=len(Y):
        raise DimException(X,Y)
    
    n = len(Y)
    result=0.0
    l = lagrange_inter_base(X,Y,x) # 计算基函数值
    for i in range(n):
        result = result + l[i]*Y[i]  #通过插值多项式计算L(x)的值
    return result

3、实例求解
课本50页 例2.2.1 设已知$f(x)=e^{-x}$的4个函数值如下表所示：

x	0	1	2	3
$e^{-x}$	1	0.3679	0.1352	0.0498

试用Lagrange公式的线性插值求$e^{-1.4}$的近似值，用二次插值求$e^{-2.1}$。

# 用Lagrange公式的线性插值求$e^{-1.4}$的近似值
print("(1) 用Lagrange公式的线性插值求$e^{-1.4}$的近似值")
X = [1,2]
Y = [0.3679,0.1353]
x = 1.4
result = lagrange_inter(X,Y,x)
print("e^{-1.4}的Lagrange插值计算结果为：",result)

返回结果如下：
(1) 用Lagrange公式的线性插值求$e^{-1.4}$的近似值
e^{-1.4}的Lagrange插值计算结果为： 0.27486000000000005

# 用二次插值求$e^{-2.1}$
print("(2) 用Lagrange公式的线性插值求$e^{-1.4}$的近似值")
X = [1,2,3]
Y = [0.3679,0.1353,0.0498]
x = 2.1
result = lagrange_inter(X,Y,x)
print("e^{-2.1}的Lagrange插值计算结果为：",result)

返回结果如下：
(2) 用Lagrange公式的线性插值求$e^{-1.4}$的近似值
e^{-2.1}的Lagrange插值计算结果为： 0.1201305

课本50页 例2.2.2 试证：由$n+1$个互异节点$x_i\in [a,b]$构成的$n+1$个插值基函数$l_i(x)(i=0,1,\cdots ,n)$具有性质：
$$\sum _{i=0}^n{l}_i(x)=1,\forall x\in [a,b]$$

# 这里采用数值方法简单验证
import numpy as np
X = [0,1,2,3]
Y = [1,0.3679,0.1353,0.0498]
x_list = np.linspace(0.1,2.9,100)
for x in x_list:
    li = lagrange_inter_base(X,Y,x)
    sum_li = sum(li)
    print(sum_li)

结果略！

5.2 Newton插值公式

5.2.1 均差及其性质

均差（也称为差商）是数值方法中的一个重要概念，它可以反映出列表函数的性质，并能对Lagrange插值公式给出新的表达形式，也就是Newton插值公式 。

• 零阶均差
  $$f[x_i]=f(x_i)$$
• 1阶均差
  $$f[x_i,x_j]=\frac{f(x_i)-f(x_j)}{x_i-x_j}$$
• 2阶均差
  $$f[x_i,x_j,x_k]=\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k}$$
• 递推地，n阶均差
  $$f[x_0,x_1,\cdots ,x_k]=\frac{f[x_0,x_1,\cdots ,x_{k-1}]-f[x_1,x_2,\cdots ,x_k]}{x_1-x_k}$$

均差具有如下特性：

• k阶均差可以表示成k+1个函数值的线性组合；
• 均差对其节点是对称的；

5.2.2 Newton插值公式

未完待续

Newton插值的Python实现

5.3 Hermite插值公式

未完待续