函数插值的python实现——拉格朗日、牛顿插值

1. 拉格朗日(Larange)插值

\qquad 拉格朗日插值，是通过构造 $n+1$ 个插值节点来构造 $n$ 阶插值多项式。

\qquad 假设 $n+1$ 个插值节点为 $x_0,x_1,\cdots ,x_n$，且满足函数插值条件：

$L_n(x_k)=f(x_k)=f_k,k=0,1,\cdots ,n$

\qquad
\qquad 拉格朗日插值的 $n$ 阶插值基函数为：

$l_k(x)=\frac{(x-x_0)\cdots (x-x_{k-1})(x-x_{k+1})\cdots (x-x_n)}{(x_k-x_0)\cdots (x_k-x_{k-1})(x_k-x_{k+1})\cdots (x_k-x_n)}$

\qquad
\qquad 显然，对于任意的插值节点 $x_i$，插值基函数满足：

$l_k(x_i)=\{\begin{array}{r}1,x_i=x_k\\ \\ 0,x_i\ne x_k\end{array}$

\qquad 因此，$n$ 阶拉格朗日插值多项式就可以写为：$L_n(x)=\sum _{k=0}^n{f}_k{l}_k(x)$

假设 $f^{(n)}(x)$ 在 $[a,b]$ 上连续，$f^{(n+1)}(x)$ 在 $(a,b)$ 上存在
假设 $n+1$ 个插值节点 $a\le x_0<x_1<\cdots <x_n\le b$
记 $n+1$ 阶多项式 ${\omega }_{n+1}(x)=(x-x_0)(x-x_1)\cdots (x-x_n)$
插值余项：
　　　　$R_n(x)=f(x)-L_n(x)=\frac{f^{(n+1)}(\xi )}{(n+1)!}{\omega }_{n+1}(x),\xi \in (a,b),\forall x\in [a,b]$

\qquad

实现代码

import numpy as np
import matplotlib.pyplot as plt

def intLagrange(x,xcp,f):
    hatf = np.zeros(len(x))
    for n in range(len(x)):
        for k in range(len(xcp)):
            xcp1 = xcp.copy()
            xcp1[k] = xcp[k] - 1
            xcp2 = xcp.copy()
            xcp2[k] = x[n] - 1
            hatf[n] += f(xcp[k])*np.prod(x[n]-xcp2)/np.prod(xcp[k]-xcp1)
    return hatf

if __name__ == '__main__':    
    plt.figure(figsize=(7, 4))    
    # parameters
    cpnum = 17      # 插值节点数
    xs = -9
    xe = 11
    x = np.arange(xs,xe,0.1)    #待插值点
    xcp = np.linspace(xs,xe,cpnum)    
    # function to be interpolated
    f = lambda x: x**2 + 30*np.sin(x)  
    plt.plot(x,f(x))    
    # Lagrange    
    hatf = intLagrange(x,xcp,f)        
    plt.plot(x,hatf)        
    # plot control point
    for i in range(len(xcp)):    
        plt.plot(xcp[i],f(xcp[i]),'ro',markerfacecolor='none')        
    plt.legend(['original','Lagrange','control point'])
    plt.title('Lagrange Interpolation')
    plt.show()

运行结果：
\qquad

拉格朗日插值法的主要缺点：一旦插值节点的数量 $n$ 发生改变，所有的插值基函数 $l_k(x),k=0,1,\cdots ,n$ 都需要重新计算。

\qquad

2. 牛顿(Larange)插值

2.1 牛顿插值多项式的基本形式

\qquad 拉格朗日插值法的主要缺点：一旦插值节点的数量 $n$ 发生改变，所有的插值基函数 $l_k(x),k=0,1,\cdots ,n$ 都需要重新计算。

\qquad 牛顿插值法可以克服这一缺点，$n$ 阶牛顿插值多项式的形式为：

$P_n(x)=a_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})$

显然，$n$ 阶牛顿插值多项式 $P_n(x)$ 只与 $n$ 个插值节点 $x_0,x_1,\cdots ,x_{n-1}$ 有关。

\qquad
\qquad 如果增加一个新的插值节点 $x_n$，那么新的 $n+1$ 阶牛顿插值多项式 $P_{n+1}(x)$ 只需要在 $P_n(x)$ 后再增加一项：

$P_{n+1}(x)=P_n(x)+$$a_{n+1}(x-x_0)(x-x_1)\cdots (x-x_{n-1})(x-x_n)$

\qquad
\qquad 对于 $n$ 阶牛顿插值多项式 $P_n(x)$，需要求出 $n$ 个插值系数 $a_0,a_1,\cdots ,a_{n-1}$，并且满足函数插值条件：

$P_n(x_i)=f(x_i),i=0,1,\cdots ,n-1$
\qquad

2.2 牛顿均差插值多项式

\qquad 将某个插值节点 $x_k$ 的函数值记为 $f_k=f(x_k)$，已知区间 $[x_0,x_{n-1}]$ 上的 $n$ 个插值节点 $x_0,x_1,\cdots ,x_{n-1}$ 的函数值 $f_0,f_1,\cdots ,f_{n-1}$。

\qquad 分析牛顿插值多项式：

$P_n(x)=a_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})$

\qquad 可以依次求出牛顿插值多项式的各个系数：

$(1)$ 当 $x=x_0$ 时，$P_n(x_0)=a_0=f_0$
$⟹a_0=f_0$

$(2)$ 当 $x=x_1$ 时，$P_n(x_1)=a_0+a_1(x_1-x_0)=f_1$
$⟹a_1=\frac{f_1-f_0}{x_1-x_0}$

$(3)$ 当 $x=x_2$ 时，$P_n(x_2)=a_0+a_1(x_2-x_0)+a_2(x_2-x_0)(x_2-x_1)=f_2$
$⟹a_2=\frac{\frac{f_2-f_0}{x_2-x_0}-\frac{f_1-f_0}{x_1-x_0}}{x_2-x_1}$
\qquad
$(4)$ 依次递推，可以求出其余系数 $a_3,\cdots ,a_{n-1}$
\qquad

(1) 均差的定义

\qquad 定义一阶均差：$f[x_0,x_k]=\frac{f(x_k)-f(x_0)}{x_k-x_0}$

\qquad 　　二阶均差：$f[x_0,x_1,x_k]=\frac{f[x_0,x_k]-f[x_0,x_1]}{x_k-x_1}$

\qquad 以及k阶均差：

$f[x_0,x_1,\cdots ,x_k]=\frac{f[x_0,\cdots ,x_{k-2},x_k]-f[x_0,x_1,\cdots ,x_{k-1}]}{x_k-x_{k-1}}$

引入均差之后，显然 $a_0=f(x_0),a_1=f[x_0,x_1],a_2=f[x_0,x_1,x_2],\cdots$

\qquad 那么 $n$ 阶牛顿插值多项式就可以写为：

$\begin{array}{rl}{N}_n(x)& =f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)\\ & +\cdots +f[x_0,x_1,\cdots ,x_n](x-x_0)(x-x_1)\cdots (x-x_{n-1})\end{array}$
\qquad

(2) 均差表

\qquad 为了简化计算，可根据均差的对称性(均差与节点的排列次序无关)，将 $k$ 阶均差的表达式写为：

$f[x_0,x_1,\cdots ,x_k]=\frac{f[x_1,x_2,\cdots ,x_k]-f[x_0,x_1,\cdots ,x_{k-1}]}{x_k-x_0}$

\qquad 根据该公式，就可以将 $n$ 阶牛顿插值多项式的插值系数 $a_0,a_1,\cdots ,a_{n-1}$ 的求解过程列成均差表：

$x_k$	$f(x_k)$	一阶均差	二阶均差	三阶均差
$x_0$	$\underline{f(x_0)}$
$x_1$	$f(x_1)$	$\underline{f[x_0,x_1]}$
$x_2$	$f(x_2)$	$f[x_1,x_2]$	$\underline{f[x_0,x_1,x_2]}$
$x_3$	$f(x_3)$	$f[x_2,x_3]$	$f[x_1,x_2,x_3]$	$\underline{f[x_0,x_1,x_2,x_3]}$
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$

在均差表中，根据 $k$ 阶均差公式，下划线标记的均差，均可由左侧的 $2$ 个(左、左上)表格项求出。

\qquad
\qquad 若记 $n$ 阶多项式 ${\omega }_n(x)=(x-x_0)(x-x_1)\cdots (x-x_{n-1}),n\ge 1$ ，那么 $n$ 阶牛顿插值多项式又可表示为：

$\begin{array}{rl}{N}_n(x)& =f(x_0)+f[x_0,x_1]{\omega }_1(x)+f[x_0,x_1,x_2]{\omega }_2(x)+\cdots +f[x_0,x_1,\cdots ,x_n]{\omega }_n(x)\end{array}$

实现代码

import numpy as np
import matplotlib.pyplot as plt

def intNewton(x,xcp,f):
    divdiff = getdivideddiff(xcp,f)
    hatf = np.ones(len(x))*f(xcp[0])
    for n in range(len(x)):
        t0 = 1
        for k in range(len(xcp)-1):
            # Newton
            t0 = (x[n]-xcp[k])*t0
            hatf[n] += divdiff[k]*t0
    return hatf

def getdivideddiff(xcp,f):
    divdiff = np.zeros(len(xcp)-1)
    t = f(xcp)
    for k in range(len(xcp)-1):
        num0 = t        
        num1 = np.concatenate((np.ones(1),num0[0:-1]))    
        k0 = k+1
        den = np.concatenate((np.ones(k0),xcp[0:-k0]))
        t = (num0-num1)/(xcp-den)
        divdiff[k] = t[k0]
    return divdiff

if __name__ == '__main__':    
    plt.figure(figsize=(7, 4))    
    # parameters
    cpnum = 17                          # 插值节点数
    xs = -9
    xe = 11
    x = np.linspace(xs,xe,100)          # 待插值点
    xcp = np.linspace(xs,xe,cpnum)      # 插值节点的函数值已知
    # function to be interpolated
    f = lambda x: x**2 + 30*np.sin(x) 
    plt.plot(x,f(x))        
    ##Newton Interpolation  
    hatf = intNewton(x,xcp,f)        
    plt.plot(x,hatf)        
    # plot control point
    for i in range(len(xcp)):    
        plt.plot(xcp[i],f(xcp[i]),'ro',markerfacecolor='none')        
    plt.legend(['original','Newton','control point'])
    plt.title('Newton Interpolation')    
    plt.show()

运行结果：
\qquad

函数插值的结果，与插值节点的选择有关。为了方便，此处使用了等距插值节点。

\qquad

3. 等距插值节点的牛顿插值

(1) 差分的定义

\qquad 考虑“插值节点之间距离相等”的情况，也就是 $x_k=x_0+kh,k=0,1,\cdots ,n$，并记 $f_k=f(x_k)$，可定义

$1)$ 一阶差分：$\{\begin{array}{r}\Delta f_k=f_{k+1}-f_k前向\\ \\ \nabla f_k=f_k-f_{k-1}后向\end{array}$

$2)$ 二阶差分：$\{\begin{array}{r}{\Delta }^2{f}_k=\Delta f_{k+1}-\Delta f_k前向\\ \\ {\nabla }^2{f}_k=\Delta f_k-\Delta f_{k-1}后向\end{array}$

\qquad 以及

$3)$ $\text{m}$阶差分：$\{\begin{array}{r}{\Delta }^m{f}_k={\Delta }^{m-1}{f}_{k+1}-{\Delta }^{m-1}{f}_k前向\\ \\ {\nabla }^m{f}_k={\Delta }^{m-1}{f}_k-{\Delta }^{m-1}{f}_{k-1}后向\end{array}$

\qquad
\qquad 从而，可得到均差与差分之间的关系：

\qquad\qquad 一阶：$f[x_k,x_{k+1}]=\frac{f_{k+1}-f_k}{x_{k+1}-x_k}=\frac{\Delta f_k}{h}$

\qquad\qquad 二阶：$f[x_k,x_{k+1},x_{k+2}]=\frac{f[x_{k+1},x_{k+2}]-f[x_k,x_{k+1}]}{x_{k+2}-x_k}=\frac{{\Delta }^2{f}_k}{2h^2}$

$⟹$$\text{m}$阶：$f[x_k,x_{k+1},\cdots ,x_{k+m}]=\frac{{\Delta }^m{f}_k}{m!h^m}$
\qquad

(2) 差分表

\qquad 等距插值节点的牛顿插值不再采用均差，而是采用差分。因此，$n$ 阶牛顿插值多项式的插值系数 $a_0,a_1,\cdots ,a_{n-1}$ 的求解过程可以列成差分表（以 $n=4$ 为例）：

$x_k$	$f_k$	$\Delta (\nabla )$	${\Delta }^2({\nabla }^2)$	${\Delta }^3({\nabla }^3)$	${\Delta }^4({\nabla }^4)$
$x_0$	$f_0$
$x_1$	$f_1$	$\underline{\Delta f_0}(\nabla f_1)$
$x_2$	$f_2$	$\Delta f_1(\nabla f_2)$	$\underline{{\Delta }^2{f}_0}({\nabla }^2{f}_2)$
$x_3$	$f_3$	$\Delta f_2(\nabla f_3)$	${\Delta }^2{f}_1({\nabla }^2{f}_3)$	$\underline{{\Delta }^3{f}_0}({\nabla }^3{f}_3)$
$x_4$	$f_4$	$\Delta f_3(\underline{\underline{\nabla f_4}})$	${\Delta }^2{f}_2(\underline{\underline{{\nabla }^2{f}_4}})$	${\Delta }^3{f}_1(\underline{\underline{{\nabla }^3{f}_4}})$	$\underline{{\Delta }^4{f}_0}(\underline{\underline{{\nabla }^4{f}_4}})$

在差分表中，单下划线标记的前向差分、双下划线标记的后向差分，均可由左侧的 $2$ 个(左、左上)表格项求出，下标正好相差一位。

\qquad

(3) 前向插值

\qquad 前向插值，是从区间 $[x_0,x_{n-1}]$ 的起点 $x_0$ 开始进行插值，任意待插值点的形式可写为 $x=x_0+th,t\ge 0$。

\qquad 牛顿插值多项式中的 $n$ 阶多项式 ${\omega }_n(x)$ 可以写为：

${\omega }_n(x)=(x-x_0)(x-x_1)\cdots (x-x_{n-1})=t(t-1)\cdots (t-n+1)h^n,t\ge 0$

\qquad
\qquad 等距插值节点时的 $n$ 阶牛顿(前向)插值多项式为：

$\begin{array}{rl}{N}_n(x)& =f(x_0)+f[x_0,x_1]{\omega }_1(x)+f[x_0,x_1,x_2]{\omega }_2(x)+\cdots \\ & +f[x_0,x_1,\cdots ,x_n]{\omega }_n(x)\\ & \\ & \\ & \\ & \Downarrow \\ & \end{array}$
$\begin{array}{rl}{N}_n(x_0+th)& =f(x_0)+\frac{\Delta f_0}{h}th+\frac{{\Delta }^2{f}_0}{2h^2}t(t-1)h^2+\cdots +\frac{{\Delta }^n{f}_0}{n!h^n}t(t-1)\cdots (t-n+1)h^n\end{array}$
$\begin{array}{rl}& =f_0+t\Delta f_0+\frac{t(t-1)}{2}{\Delta }^2{f}_0+\cdots +\frac{t(t-1)\cdots (t-n+1)}{n!}{\Delta }^n{f}_0\end{array}$

写成该表达式的目的，是为了更方便地使用“差分表”来进行实现。

\qquad

(4) 后向插值

\qquad 后向插值，正好与前向插值相反，是从区间 $[x_0,x_{n-1}]$ 的终点 $x_{n-1}$ 开始进行插值，任意待插值点的形式可写为 $x=x_{n-1}+th,t\le 0$。为了能够利用前向插值公式，可将插值节点按照 $x_{n-1},x_{n-2},\cdots ,x_1,x_0$ 的顺序进行倒序排列。

\qquad 倒序排列后，$n$ 阶牛顿插值多项式为：

$\begin{array}{rl}{N}_n(x)& =f(x_{n-1})+f[x_{n-1},x_{n-2}](x-x_{n-1})\\ & +f[x_{n-1},x_{n-2},x_{n-3}](x-x_{n-1})(x-x_{n-2})+\cdots \\ & +f[x_0,x_1,\cdots ,x_n](x-x_{n-1})(x-x_{n-2})\cdots (x-x_0)\end{array}$

\qquad 同样地，可以记 $n$ 阶多项式 ${\omega }_n(x)$ 为：

${\omega }_n(x)=(x-x_{n-1})(x-x_{n-2})\cdots (x-x_0)=t(t+1)\cdots (t+n-1)h^n$

\qquad 以及 $\text{m}$阶差分：$f[x_k,x_{k-1},\cdots ,x_{k-m}]=\frac{{\nabla }^m{f}_k}{m!h^m}$

\qquad
\qquad 因此，等距插值节点时的 $n$ 阶牛顿(后向)插值多项式为：

$\begin{array}{rl}{N}_n(x)& =f(x_{n-1})+f[x_{n-1},x_{n-2}]{\omega }_1(x)+f[x_{n-1},x_{n-2},x_{n-3}]{\omega }_2(x)+\cdots \\ & +f[x_{n-1},x_{n-2},\cdots ,x_0]{\omega }_n(x)\\ & \\ & \\ & \\ & \Downarrow \\ & \end{array}$
$\begin{array}{rl}{N}_n(x_0+th)& =f(x_{n-1})+\frac{\nabla f_{n-1}}{h}th+\frac{{\nabla }^2{f}_{n-1}}{2h^2}t(t-1)h^2+\cdots +\frac{{\nabla }^n{f}_{n-1}}{n!h^n}t(t-1)\cdots (t-n+1)h^n\end{array}$
$\begin{array}{rl}& =f_{n-1}+t\nabla f_{n-1}+\frac{t(t+1)}{2}{\nabla }^2{f}_{n-1}+\cdots +\frac{t(t+1)\cdots (t+n-1)}{n!}{\nabla }^n{f}_{n-1}\end{array}$

写成该表达式的目的，是为了更方便地使用“差分表”来进行实现。

\qquad

实现代码

import numpy as np
import matplotlib.pyplot as plt

# 牛顿法：等距节点插值（前插公式）
def intNewtonforw_e(x,xcp,f,h,num):     
    diff_forw, diff_back = getdiff(xcp,f)
    hatf = np.zeros(len(x))
    for i in range(len(x)):
        if i%num == 0:
            hatf[i] = f(x[i])
            continue
        t = i/num
        hatf[i] = hatf[0]
        tmp = 1
        prodn = 1
        for n in range(1,len(xcp)-1):#插值多项式阶数为n=len(xcp)-1阶，可以小于n
            tmp = tmp*(t-n+1)
            prodn = prodn * n
            hatf[i] += tmp * diff_forw[n]/prodn                   
    return hatf

# 牛顿法：等距节点插值（后插公式）
def intNewtonback_e(x,xcp,f,h,num):     
    diff_forw, diff_back = getdiff(xcp,f)
    xn = len(x)-1
    hatf = np.zeros(len(x))
    for i in range(len(x)):
        if i%num == 0:
            hatf[xn-i] = f(x[xn-i])
            continue
        t = - i/num
        hatf[xn-i] = hatf[xn]
        tmp = 1
        prodn = 1
        for n in range(1,len(xcp)-1):#插值多项式阶数为nlen(xcp)-1阶，可以小于n
            tmp = tmp*(t+n-1)
            prodn = prodn * n
            hatf[xn-i] += tmp * diff_back[n]/prodn                     
    return hatf
    
# 获取各阶前向差分和后向差分    
def getdiff(xcp,f):
    # np.set_printoptions(precision=5)
    diff_forw = np.zeros(len(xcp))
    diff_back = np.zeros(len(xcp))
    t0 = f(xcp)
    diff_forw[0] = t0[0]
    diff_back[0] = t0[-1]    
    for n in range(1,len(xcp)):
        fn0 = t0        
        fn1 = np.concatenate((np.ones(1),fn0[0:-1]))    
        t1 = fn0 - fn1
        diff_forw[n] = t1[1]
        diff_back[n] = t1[-1]
        t0 = t1[1:]
    return diff_forw, diff_back        

if __name__ == '__main__':    
    plt.figure(figsize=(7, 4))    
    # parameters
    xs = -13
    xe = 13
    h = 1                           # 相邻等距节点的间隔
    num = 10                        # 两个等距节点之间再等距划分10个点
    t = h/num                       # 相邻节点之间的待插值点间隔
    x = np.arange(xs,xe+t,t)        # 等距节点之间再等分后的待插值点
    xcp = np.arange(xs,xe+h,h)      # 等距节点
    # function to be interpolated
    f = lambda x: x**2 + 30*np.sin(x)     
    plt.plot(x,f(x))        
    ##Newton Interpolation  
    hatf = intNewtonforw_e(x,xcp,f,h,num) #前向插值
    plt.plot(x,hatf) 
    hatf = intNewtonback_e(x,xcp,f,h,num) #后向插值   
    plt.plot(x,hatf)        
    # plot control point
    for i in range(len(xcp)):    
        plt.plot(xcp[i],f(xcp[i]),'ro',markerfacecolor='none')        
    plt.legend(['original','Newton forward','Newton backward','control point'])
    plt.title('Newton Interpolation')
    plt.show()

运行结果：

$(1)$ 前插公式
\qquad

$(2)$ 后插公式
\qquad

$(3)$ 将前插公式和后插公式的结果显示在同一张图
\qquad
\qquad
参考：
基于径向基函数(RBF)的函数插值
\qquad
\qquad
\qquad
\qquad