Lagrange插值

# Lagrange插值

标签（空格分隔）： 数值方法 插值

---

## 插值

在实际问题中，往往只能观察到某物理量在有限个点处的值。如何从这有限个点值将物理量的函数恢复出来，进而研究该物理量的特性呢？

另外，在计算机中，只能处理有限、离散的信息。它可以通过保存一些系数来表达特定的函数，如：通过存贮$a_0,a_1,\cdots,a_n$来表示多项式$a_0+a_1x+\cdots+a_nx^n$，

或者存贮$a,b,c$来表示函数$a\sin(bx+c)$。但对于一个任意的函数$f(x)$，计算机就没办法来表示。一种可选的方法是存贮函数$f(x)$在有限个点处的函数值，然后在需要的时候，尽可能的将函数$f(x)$恢复出来。

它们的共同的**数学模型**是：已知函数$f(x)$在若干个点处的函数值，如何恢复这个原来的函数$f(x)$？

1. 通过有限个点的函数有无穷多个。因此，想要得到原来的函数是几乎不可能的；

2. 只能找一个原来函数的近似函数。因此，可以找一个性质“好”的，形式上熟悉的近似函数；

$\{x_i\}_{i=0}^n$是$n+1$个互不相同的点，函数$f(x)$在这些节点上的函数值为$f(x_i)$，$\Phi$是已知的函数类(函数空间)。若$\Phi$中存在一个函数$g(x)$，使得

$$

g(x_i)=f(x_i), i=0,1,\cdots,n

$$

则称$g(x)$是$f(x)$的一个*插值函数*。$f(x)$称为*原函数*。$\{x_i\}$称为*插值节点*。

找一个性质“好”的近似函数，可以通过使用熟悉的函数类来得到。如：

1. 多项式 $P_n=span\{1,x,x^2,\cdots,x^n\}$；

2. 三角函数 $S_n=span\{1,\cos x,\sin x,\cos 2x, \sin 2x,\cdots,\cos nx,\sin nx\}$

## 插值函数的存在唯一性

怎么得到$g(x)$？

设$\Phi=span\{\phi_0(x),\cdots,\phi_m(x)\}$是一个$m+1$维的线性空间，则

$$

g(x)=a_0\phi_0(x)+\cdots+a_m\phi_m(x)

$$

想办法求出$a_0,\cdots,a_m$即找到了$g(x)$。由插值的要求，$g(x)$满足

$$

g(x_i)=a_0\phi_0(x_i)+\cdots+a_m\phi_m(x_i)=f(x_i), i=0,1,\cdots,n

$$

写成矩阵形式，有

$$

\begin{pmatrix}

\phi_0(x_0) & \phi_1(x_0) & \cdots & \phi_m(x_0) \\

\phi_0(x_1) & \phi_1(x_1) & \cdots & \phi_m(x_1) \\

\vdots & \vdots & & \vdots \\

\phi_0(x_n) & \phi_1(x_n) & \cdots & \phi_m(x_n) \\

\end{pmatrix}

\begin{pmatrix}

a_0 \\ a_1 \\ \vdots \\ a_m

\end{pmatrix}

=\begin{pmatrix}

f(x_0) \\ f(x_1) \\ \vdots \\ f(x_n)

\end{pmatrix}

$$

这是一个关于系数$a_0, \cdots, a_m$的线性方程组。 利用线性代数相关知识，可以分析出这个方程的解结构。特别地，有

**定理**

若$m=n$，且

$$

\left|\begin{matrix}

\phi_0(x_0) & \phi_1(x_0) & \cdots & \phi_n(x_0) \\

\phi_0(x_1) & \phi_1(x_1) & \cdots & \phi_n(x_1) \\

\vdots & \vdots & & \vdots \\

\phi_0(x_n) & \phi_1(x_n) & \cdots & \phi_n(x_n) \\

\end{matrix}\right|\neq 0

$$

则插值函数$g(x)$存在唯一。

## 插值多项式

在已知的函数形态中，多项式是最早接触到、形式简单的一类函数。它特别适合计算机使用

1. 计算一个多项式只需要四则运算，方便编程

2. 一个多项式的积分与导数仍然是一个多项式，且可以方便的得到任意阶导数和任意次的积分

多项式还有一些特性，如：

1. $f(x)$是一个$n$次多项式，若$f(a)=0$，则有$f(x)=g(x)(x-a)$，且$g(x)$是一个$n-1$次多项式

2. 若$f(a)=0$，且$f'(a)=0$，则有$f(x)=g(x)(x-a)^2$，且$g(x)$是一个$n-2$次多项式

由于多项式的诸多优秀的特性，在以后的学习中，一直使用多项式来近似。

回到插值过程，若节点$\{x_i\}_{i=0}^n$互不相同，在$n$次多项式空间$P_n=span\{1,x, x^2,\cdots, x_n\}$中找插值函数。由线性代数知识，有

$$

\left|\begin{matrix}

1 & x_0 & \cdots & x_0^n \\

1 & x_1 & \cdots & x_1^n \\

\vdots & \vdots & & \vdots \\

1 & x_n & \cdots & x_n^n

\end{matrix}\right|=\prod_{i,j=0, \\ j>i}^n(x_j-x_i)\neq 0

$$

这是一个Vandoormod行列式。因此，有结论

**定理**

$n+1$个互不相同节点上的$n$次插值多项式存在且唯一。

## 待定系数法得到插值多项式

前面已经给出了求解插值多项式的一种方法--待定系数。

设插值多项式是$g(x)=a_0+a_1x+\cdots+a_n x^n$，由插值条件，可以得到线性方程组

$$

\begin{cases}

a_0+a_1 x_0+\cdots+a_n x_0^n & = f(x_0) \\

a_0+a_1 x_1+\cdots+a_n x_1^n & = f(x_1) \\

\cdots \\

a_0+a_1 x_n+\cdots+a_n x_n^n & = f(x_n) \\

\end{cases}

$$

解出系数$a_0, \cdots, a_n$即可。

**例**

求过点$(1, 3)$, $(2, 2)$, $(3, 4)$的插值多项式。

**解**

总共3个点，因此使用2次多项式。设$g(x)=a_0+a_1 x+a_2 x^2$，则有

$$

\begin{cases}

a_0+ a_1 + a_2  &=3 \\

a_0+2a_1 + 4 a_1&=2 \\

a_0+3a_1 + 9 a_2&=4 \\

\end{cases}

$$

**方法的不足**：需要解一个线性方程组，而且Vandoormod矩阵是一个病态矩阵。当$n$比较大($>10$)时，数值上的小扰动，会带来明显的计算误差。

## Lagrange插值

设插值多项式是

$$

g(x)=f(x_0)l_0(x)+f(x_1)l_1(x)+\cdots+f(x_n)l_n(x)

$$

其中$l_i(x), i=0,\cdots,n$是次数不超过$n$次的多项式。下面，想办法救出$l_i(x)$。

由

$$

f(x_0)=g(x_0)=f(x_0)l_0(x_0)+f(x_1)l_1(x_0)+\cdots+f(x_n)l_n(x_0)

$$

因此，可以假定成立

$$

l_0(x_0)=1, l_i(x_0)=0, i=1,2,\cdots,n

$$

类似地，可以得到

$$

l_0(x_1)=0, l_1(x_1)=1, l_i(x_1)=0, i=2,3,\cdots,n

$$

对$i=0,1,\cdots,n$，都有

$$

l_i(x_i)=1, l_j(x_i)=0, j=0,1,\cdots,j-1,j+1,\cdots,n

$$

可以得到如下的表

-  | $l_0(x)$ | $l_1(x)$ | $\cdots$ | $l_n(x)$

- | :-: | :-: | :-: | -: |

| $x_0$  | 1 | 0 | $\cdots$ | 0

| $x_1$ | 0 | 1 | $\cdots$ | 0

| $\vdots$ | $\vdots$

|$x_n$ | 0 | 0 | $\cdots$ |1

可以看出，每个函数$l_i(x)$满足

$$

l_i(x_j)=\delta_{ij}=\begin{cases} 1, & i=j \\ 0, & i\neq j \end{cases}

$$

$l_i(x)$有零点$x_0$, $x_1$, $\cdots$, $x_{i-1}$, $x_{i+1}$, $\cdots$, $x_n$， 因此有

$$

l_i(x)=a_i(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)

$$

其中$a_i$是某个实数。又$l_i(x_i)=1$，即

$$

a_i(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)=1

$$

解得

$$

a_i=\frac1{(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)}

$$

得到

$$

l_i(x)=\frac{(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)}

$$

称$l_i(x)$为*Lagrange基函数*。插值函数

$$

g(x)=f(x_0)l_0(x)+f(x_1)l_1(x)+\cdots+f(x_n)l_n(x)

$$

为*Lagrange插值*，记为$L_n(x)$。

**例**

求过点$(1, 3)$, $(2, 2)$, $(3, 4)$的插值多项式。

**解**

插值节点是$1,2,3$，因此Lagrange插值基函数是

$$

\begin{aligned}

l_0(x)=&\frac{(x-2)(x-3)}{(1-2)(1-3)} \\

l_1(x)=&\frac{(x-1)(x-3)}{(2-1)(2-3)} \\

l_2(x)=&\frac{(x-1)(x-2)}{(3-1)(3-2)} \\

\end{aligned}

$$

因此，插值函数为

$$

g(x)=3\frac{(x-2)(x-3)}{(1-2)(1-3)}+2\frac{(x-1)(x-3)}{(2-1)(2-3)}

+4\frac{(x-1)(x-2)}{(3-1)(3-2)}

$$

可以看到$L_n(x)$是一个次数不超过$n$次的多项式，由插值多项式的存在唯一性知，$L_n(x)$与待定系数法得到的多项式是一样的，而它避开了求解一个病态的线性方程组。

## Lagrange插值的伪代码

``` c

fx=0.0;

for i=0,n:

  li=1.0;

  for j=0,i-1:

    li=li*(t-x[j])/(x[i]-x[j]);

  end for

  for j=i+1,n:

    li=li*(t-x[j])(x[i]-x[j])

  end for

  fx = fx+y[i]*li;

end for

// fx 即为 插值函数 在 t 处的值

```

## 误差分析

插值多项式当然跟原函数有误差，下面来分析一下这个误差。记

$$

R_n(x)=f(x)-L_n(x)

$$

则，由插值的定义知

$$

R_n(x_i)=f(x_i)-L_n(x_i)=0, i=0,1,\cdots,n

$$

因此，可以设

$$

R_n(x)=K_n(x)(x-x_0)(x-x_1)\cdots(x-x_n)

$$

其中$K_n(x)$是一个未知的函数。

下面，想办法把$K_n(x)$表达出来。对$f(x)$定义域内任意的数$a$，构造辅助函数

$$

\phi(t)=f(t)-L_n(t)-K_n(a)(t-x_0)\cdots(t-x_n)

$$

为$t$的函数，则有

$$

\phi(x_i)=f(x_i)-L_n(x_i)=0, i=0,1,\cdots,n

$$

及

$$

\phi(a)=f(a)-L_n(a)-R_n(a)=0

$$

也就是说，点$x_i$, $i=0,\cdots,n$及$a$均为函数$\phi(t)$的零点。设所有这些点在区间$[b,c]$内，则

* $\phi(t)$具有$n+2$个互不相同的零点，若$\phi(t)$可导，则由Rolle定理，至少存在$n+1$个互不相同的点是$\phi'(t)$的零点，这些零点在区间$(b,c)$内。

* 同样的，若$\phi'(t)$可导，则至少存在$(b,c)$内的$n$个互不相同的$\phi''(t)$的零点。

* 因此，若$\phi(t)$具有$n+1$阶导数，则至少存在一个点$\xi\in(b,c)$满足$\phi^{(n+1)}(\xi)=0$。

* 从$\phi(t)$的表达式中可以看到，$\phi(t)$的光滑性与$f(t)$的光滑性相关。若$f(x)$足够光滑（比如：有到$m$阶的连续导数），则$\phi(t)$也具有$m$阶连续导数。

因此，若$f(x)$具有$n+1$阶导数，则存在$\xi$满足

$$

0=\phi^{(n+1)}(\xi)=f^{(n+1)}(\xi)-K_n(a)(n+1)!

$$

即有

$$

K_n(a)=\frac{f^{(n+1)}(\xi)}{(n+1)!}

$$

由$a$的任意性，可以得到误差表达式

$$

R_n(x)=f(x)-L_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^n(x-x_i)

$$

**例**

已知$\sin\frac{\pi}6=\frac12$, $\sin\frac{\pi}4=\frac{\sqrt 2}2$, $\sin\frac{\pi}3=\frac{\sqrt 3}2$，用1次和2次Lagrange插值近似$\sin50^\circ$，并估计误差

**解**

求解位置$x=50^\circ=\frac{5\pi}{18}$。记$x_0=\frac{\pi}6$, $x_1=\frac{\pi}4$, $x_2=\frac{\pi}3$。

1. 先用$x_0$, $x_1$做线性插值，近似计算$x$处的值

$$

L_1(x)=\frac12\frac{x-\frac{\pi}4}{\frac{\pi}6-\frac{\pi}4}

+\frac{\sqrt2}2\frac{x-\frac{\pi}6}{\frac{\pi}4-\frac{\pi}6}

$$

近似得到$\sin50^\circ\approx L_1(\frac{5\pi}{18})\dot=0.77614$。根据误差表达式

$$

R_1(x)=\frac{\sin^{(2)}(\xi)}{2!}(x-\frac{\pi}6)(x-\frac{\pi}4),\quad  \xi\in(\frac{\pi}6,\frac{\pi}3)

$$

又

$$

\sin^{(2)}(\xi)=-\sin(\xi)\in(-\frac{\sqrt3}2,-\frac12), \quad \xi\in(\frac{\pi}6,\frac{\pi}3)

$$

得到

$$

R_1(\frac{5\pi}{18})=-\sin(\xi)\frac12\frac{\pi}{18}\frac{4\pi}{18}

\in(-0.01319,-0.00762)

$$

比较真解$\sin50^\circ=0.7660444\cdots$，得到误差约为$-0.01001$，落在估计的区间内。

2. 类似地，用$x_1$, $x_2$做插值，得到估计值和误差$\tilde R_1$

$$

\sin50^\circ\approx0.76008, \quad \tilde R_1\in(0.00538, 0.00660)

$$

实际误差约为$0.00596$

3. 用$x_0$, $x_1$, $x_2$三个点构造二次插值多项式

$$

L_2(x)=\frac12\frac{(x-\frac{\pi}4)(x-\frac{\pi}3)}{(\frac{\pi}6-\frac{\pi}4)(\frac{\pi}6-\frac{\pi}3)}

+\frac{\sqrt2}2\frac{(x-\frac{\pi}6)(x-\frac{\pi}3)}{(\frac{\pi}4-\frac{\pi}6)(\frac{\pi}4-\frac{\pi}3)}

+\frac{\sqrt3}2\frac{(x-\frac{\pi}6)(x-\frac{\pi}4)}{(\frac{\pi}3-\frac{\pi}6)(\frac{\pi}3-\frac{\pi}4)}

$$

及误差

$$

R_2(x)=\frac{\sin^{(3)}(\xi)}{3!}(x-\frac{\pi}6)(x-\frac{\pi}4)(x-\frac{\pi}3), \xi\in(\frac{\pi}6,\frac{\pi}3)

$$

得到

$$

L_2(\frac{5\pi}{18})\approx0.76543, \quad R_2(\frac{5\pi}{18})\in(0.00044,0.00077)

$$

实际误差约为$0.00061$

* 在插值中，所求点$x$落在节点$\{x_i\}_{i=0}^n$所组成的区间外部的插值，又称为*外插*。从上面的例子可以看到，通常内插的误差要小于外插的误差。

* 上面的例子还可以看到，高阶插值的误差要小于低阶插值的误差。

## 事后误差估计

在误差表达式中，$f^{(n+1)}(\xi)$在实际应用中，是不可计算的。

1. 若$f(x)$已知，则可以通过$\xi$的范围来得到$f^{(n+1)}(\xi)$的取值范围，进而给出误差的取值范围来。

2. 若$f(x)$也未知，则没办法得到误差的范围。

在实际应用中，**事后误差估计**方法可以对误差进行估计。

设节点$\{x_0,x_1,\cdots,x_n\}$上的插值多项式是$L_n(x)$，则有误差

$$

f(x)-L_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}(x-x_0)\cdots(x-x_n)

$$

替换一个节点，节点组$\{x_1,x_2,\cdots,x_{n+1}\}$上的插值多项式是$\tilde L_n(x)$，误差为

$$

f(x)-\tilde L_n(x)=\frac{f^{(n+1)}(\eta)}{(n+1)!}(x-x_1)\cdots(x-x_{n+1})

$$

这里，做一个假设$f^{(n+1)}(\eta)\approx f^{(n+1)}(\xi)$，则有

$$

\frac{f(x)-L_n(x)}{f(x)-\tilde L_n(x)}\approx\frac{x-x_0}{x-x_{n+1}}

$$

得到

$$

f(x)\approx\frac{x-x_{n+1}}{x_0-x_{n+1}}L_n(x)+\frac{x-x_0}{x_{n+1}-x_0}\tilde L_n(x)

$$

进而，$L_n(x)$的误差可以估计为

$$

f(x)-L_n(x)\approx\frac{x-x_0}{x_{n+1}-x_0}\left(\tilde L_n(x)-L_n(x)\right)

$$

这样，得到了**事后误差估计**。

**注意**： 推导过程中，假设$f^{(n+1)}(\eta)\approx f^{(n+1)}(\xi)$ 是没有数学依据的， 因此得到的误差估计即可能比真实误差大，也可能比真实误差小。从数学上说，这个数值没有意义。

但实际使用中，效果挺好。