插值方法以及Runge现象

在数值分析中，插值方法是基础且重要的。本文将介绍Lagrange插值公式与Newton插值公式。此外，针对Runge现象，本文给出了稍稍详细的讨论。

一、Lagrange插值公式

假设函数 \(y=f(x)\) 在取定的\(n+1\)个互异的基点 \(x_0,x_1,\cdots,x_n\) 处的值已知分别为 \(y_0=f(x_0),y_1=f(x_1),\cdots,y_n=f(x_n)\)，现在要寻找多项式 \(p(x)\) 使得 $$p(x_k)=f(x_k), \quad k=0,1,\cdots,n$$
记

\[l_i(x)=\prod_{j=0, j \neq i}^n \frac{x-x_j}{x_i-x_j},\quad i=0,1,\cdots,n \]

令

\[p_n(x)=\sum_{i=0}^nf(x_i)l_i(x) \]

这就是满足条件的多项式。

例 利用等距基点计算函数 \(f(x)=\frac{\sin{e^{-x^2}}}{\sqrt{x+2}}\) 在区间 \([-1,1]\) 上的值。

import math
def fun1(t):
    return math.sin(math.exp(-t**2))/math.sqrt(t+2)

def LagIn(X,P,x):
    Pol=0
    length=len(P)
    count1=0
    while count1 <= length-1:
        lpoly=P[count1]
        count2=0
        while count2 <= length-1:
            if count1 != count2:
                lpoly=lpoly*(x-X[count2])/(X[count1]-X[count2])
            else:
                lpoly=lpoly
            count2=count2+1
        Pol=Pol+lpoly
        count1=count1+1
    return Pol    

n=int(input("输入基点个数："))
A=[]
K=[]
Result=[]
count=0
while count <= n:
    A.append(fun1(-1+count*2/n))
    K.append(-1+count*2/n)
    count=count+1
count=0
print(LagIn(K,A,0.1))

考虑误差问题。记

\[r_n(x)=f(x)-p_n(x) \]

若\(f(x)\)在包含\(n+1\)个插值基点\(x_0,x_1,\cdots,x_n\)的区间\([a,b]\)上有\(n\)阶连续导数，且在\((a,b)\)内存在\(n+1\)阶有界导数，那么对\([a,b]\)上的每一点\(x\)必存在一点\(\eta \in (a,b)\)使得

\[r_n(x)=\frac{f^{(n+1)}(\eta)}{(n+1)!}\prod_{i=0}^n(x-x_i) \]

二、均差与Newton插值公式

仍然使用上一节的记号。考虑形如

\[N_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_n(x_i)=f(x_i),\quad i=0,1,\cdots,n \]

由线性代数的知识容易知道，\(a_0,a_1,\cdots,a_n\)是唯一确定的。考虑这种形式的多项式的好处是，当需要增加基点时，可以减少计算量。

下面引入均差的概念来确定该插值多项式的系数。
\(f(x)\)关于基点\(x_0,x_1,\cdots,x_n\)的n阶均差定义为

\[f[x_0,x_1,\cdots,x_n]=\frac{f[x_1,x_2,\cdots,x_n]-f[x_0,x_1,\cdots,x_{n-1}]}{x_n-x_0} \]

规定

\[f[x_i]=f(x_i) \]

是\(f(x)\)关于基点\(x_i\)的零阶均差。

利用数学归纳法可以证明：n阶均差可以写成

\[f[x_0,x_1,\cdots,x_n]=\sum_{i=0}^n\frac{f(x_i)}{\prod_{j=0,j\neq i}^n(x_i-x_j)} \]

与利用Lagrange插值法得到的多项式比较可知，\(N_n(x)\)的系数

\[a_k=f[x_0,x_1,\cdots,x_k],\quad k=0,1,\cdots,n \]

三、Runge现象

考虑函数

\[f(x)=\frac{1}{1+x^2},\quad x\in [-5,5] \]

对该函数在等距基点上插值，分别取插值基点个数为5，10，15，20，30，作出误差图像，如下所示：

import math
import matplotlib.pyplot as mpl
def fun1(t):
    return (1+t**2)**(-1)

def LagIn(X,P,x):
    Pol=0
    length=len(P)
    count1=0
    while count1 <= length-1:
        lpoly=P[count1]
        count2=0
        while count2 <= length-1:
            if count1 != count2:
                lpoly=lpoly*(x-X[count2])/(X[count1]-X[count2])
            else:
                lpoly=lpoly
            count2=count2+1
        Pol=Pol+lpoly
        count1=count1+1
    return Pol    

n=int(input("输入基点个数："))
A=[]
K=[]
Result=[]
count=0
while count <= n:
    A.append(fun1(-5+count*10/n))
    K.append(-5+count*10/n)
    count=count+1
count=0
B=[]
F=[]
while count <= 10000:
    B.append(-5.01+count*10.02/10000)
    F.append((fun1(B[count])-LagIn(K,A,B[count])))
    count=count+1

mpl.plot(B,F)
mpl.show()

可以看到，插值多项式在端点附近有较大的误差。事实上可以证明

\[\lim_{n \to \infty}\max_{x \in [-5,5]}|f(x)-p_n(x)|= \infty \]

这就是所谓的Runge现象。下面将进行仔细讨论。

在区间\([a,b]\)上，记插值基点为\(\{x_j^{(n)}\},\quad j=0,1,\cdots,n\)，插值多项式为\(p_n(x)\)。首先估计误差项\(f(x)-p_n(x)\)。
由Newton插值法，在区间\([a,b]\)上取与插值基点相异的一个基点 \(x=t\)，得到

\[p_{n+1}(x)=p_n(x)+f[x_0,x_1,\cdots,x_n,t](x-x_0)(x-x_1)\cdots(x-x_n) \]

由于\(p_{n+1}(t)=f(t)\)，因此

\[f(x)-p_n(x)=f[x_0,x_1,\cdots,x_n,x]\prod_{j=0}^n(x-x_{j}^{(n)}) \]

下面记

\[\omega_n(x)=\prod_{j=0}^n(x-x_{j}^{(n)}) \]

下面在复平面上考虑问题，插值基点视作实轴上的点，仍保留原来的记号。

由留数定理

\[f[x_0,x_1,\cdots,x_n]=\sum_{i=0}^n\frac{f(x_i)}{\prod_{j=0,j\neq i}^n(x_i-x_j)}=\sum_{i=0}^n\mathrm{Res}(\mathrm{g},x_i)=\frac{1}{2\pi i}\int_Cg(z)\mathrm{d}z \]

其中

\[g(z)=\frac{f(z)}{\prod_{i=0}^n(z-x_i)} \]

C是曲线，使得插值基点在它内部，\(f\)在其内部解析。类似可以求得\(f[z,x_0,x_1,\cdots,x_n]\)的表达式。从而得到

\[f(z)-p_n(z)=\frac{1}{2\pi i}\int_C\frac{\omega_n(z)}{\omega_n(\eta)}\frac{f(\eta)}{\eta-z}\mathrm{d}\eta \]

为了进一步估计，需要先证明以下引理：

引理 假设\(\{x_j^{(n)}\}\)是区间\([a,b]\)上的等距基点，定义

\[\sigma_n(z)=|\omega_n(z)|^{\frac{1}{n+1}} \]

则

\[\lim_{n\to \infty}\sigma_n(z)=\sigma(z)=\mathrm{exp}\{\frac{1}{b-a}\int_a^b \mathrm{ln}|z-s|\mathrm{d}s\} \]

证明：

\[\mathrm{ln}\sigma_n(z)=\frac{1}{n+1}\sum_{i=0}^n \mathrm{ln}|z-a-\frac{b-a}{n}i| \]

由黎曼积分的定义知

\[\lim_{n\to \infty}\mathrm{ln}\sigma_n(z)=\frac{1}{b-a}\int_a^b\mathrm{ln}|z-s|\mathrm{d}s \]

引理成立。

考虑曲线

\[C(\rho)=\{z \in \mathbb{C}|\sigma(z)=\rho \} \]

命题 \(\{x_j^{(n)}\}\)是包含在曲线\(C(\rho)\)中的等距插值基点，若\(f\)在\(C(\rho)\)中解析，则\(p_n\)在\(C(\rho'),\rho' < \rho\)上一致收敛于\(f\)。
证明：假设\(z\)在曲线\(C(\rho')\)上，显然\(C(\rho)\)和\(C(\rho')\)是不交的，因此有以下估计：

\[|f(z)-p_n(z)|\leq\frac{1}{2\pi i}\int_{C(\rho)}|\frac{\omega_n(z)}{\omega_n{(\eta)}}||\frac{f(\eta)}{\eta-z}|\mathrm{d}\eta \]

由于\(|f(\eta)|\)和\(|\eta-z|\)在\(C(\rho)\)上分别由上界和下界，因此

\[|f(z)-p_n(z)|\leq A \max_{\eta \in C(\rho)}|\frac{\omega_n(z)}{\omega_n(\eta)}| \]

考虑引理中的结果，由\(\lim_{n \to \infty}\sigma_n(z)=\rho'\)及\(\lim_{n \to \infty}\sigma_n(\eta)=\rho\)，对\(\epsilon < \frac{|p-p'|}{3}\)能够取到\(n\)使得

\[|\frac{\omega_n(z)}{\omega_n(\eta)}|^\frac{1}{n+1}=\frac{\sigma_n(z)}{\sigma_n(\eta)}\leq \frac{\rho'+\epsilon}{\rho-\epsilon}<1 \]

从而可以证明一致收敛。

考虑\(f\)在\(C(\rho)\)中并不解析的情况。仅仅讨论一种特殊情况：若只有一个一阶极点\(z^*\)，此时

\[f(z)-p_n(z)=\frac{1}{2\pi i}\int_{C(\rho)}\frac{\omega_n(z)}{\omega_n(\eta)}\frac{f(\eta)}{\eta-z}\mathrm{d}\eta+\frac{1}{2\pi i}\int_{C^*}\frac{\omega_n(z)}{\omega_n(\eta)}\frac{f(\eta)}{\eta-z}\mathrm{d}\eta \]

其中\(C^*=C(\rho^*)\)是围绕极点的一条路径。由上面的讨论知，第一项趋于零。利用留数定理，第二项可以写成

\[B\frac{\omega_n(z)}{\omega_n(z^*)} \]

其中

\[B=\lim_{\eta \to z^*}\frac{f(\eta)}{(\eta-z^*)(\eta-z)} \]

是有界的。

\[|\frac{\omega_n(z)}{\omega_n(z^*)}|^\frac{1}{n+1}=\frac{\sigma_n(z)}{\sigma_n(z^*)}\geq\frac{\rho'-\epsilon}{\rho^*+\epsilon}>1 \]

此时\(f(z)-p_n(z)\)并不收敛。

上面的讨论并没有给出Runge现象的充要条件。当插值区间过大时，可能会出现Runge现象；缩小插值区间，可以避免Runge现象，因此可以考虑分段插值的方法。此外也可以采用Chebyshev基点进行插值，这就不在本文的讨论范围内了。

---

参考
[1]Epperson J. On the Runge Example[J]. The American Mathematical Monthly, 1987, 94(4):329-341.
[2]林成森. 数值计算方法[M]. 科学出版社, 2005.
[3]作者本人的实验报告.