11. SciPy拉格朗日插值

拉格朗日插值法是以法国十八世纪数学家约瑟夫·路易斯·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律，而不少函数都只能通过实验和观测来了解。如对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日（插值）多项式。

11.1 拉格朗日多项式

一般下，若已知函数$\;y = f(x)$在[a, b]上有定以，且有$\;n+1$个互异的值$\;x_0,x_1,\cdots,x_n$处的函数值为$\;y_0,y_1,\cdots,y_n$,则可以考虑构造一个过这$\;n+1 $个点的、次数不超过$n$的多项式$\;y = P_n(x)$,使其满足： $$ P_n(x_k) = y_k, k \in[0, n] $$ 函数$\;y = P_n(x)$称为$\;f(x)$的插值函数，$x_0,x_1,\cdots,x_n$为插值点。设集合$D_n= ${$0,1,\cdots,n$}，对于任意$\;k \in D_n$都有$\;P_k(x)$ ,$\;B_k = ${$\;i\;\vert\; i \neq k\;,\;i \in D_n$}得到拉格朗日基本多项式(插值基函数)公式： $$ P_k(x) = \prod_{i \in B_k} \frac{x-x_i}{x_k - x_i} $$ 应用拉格朗日插值基函数得到拉格朗日插值多项式： $$ L_n(x)=\sum_{j = 0}^{n}y_jP_j(x) $$ 例如当$\;n=4$时，$D_n={0,1,2,3}$， $$ f(x) =\sum_{j = 0}^{3}y_jP_j(x)= y_0P_0(x) +y_1P_1(x) +y_2P_2(x) +y_3P_3(x) $$ 其中，$P_k(x),\;k\in D_n$为： $$ P_0(x) = \prod_{i \in B_k} \frac{x-x_i}{x_0 - x_i} = \frac{x - x_1}{x_0-x_1}\times \frac{x - x_2}{x_0-x_2}\times \frac{x - x_3}{x_0-x_3} $$ $$ P_1(x) = \prod_{i \in B_k} \frac{x-x_i}{x_1 - x_i} = \frac{x - x_0}{x_1-x_0}\times \frac{x - x_2}{x_1-x_2}\times \frac{x - x_3}{x_1-x_3} $$ $$ P_2(x) = \prod_{i \in B_k} \frac{x-x_i}{x_2 - x_i} = \frac{x - x_0}{x_2-x_0}\times \frac{x - x_1}{x_2-x_1}\times \frac{x - x_3}{x_2-x_3} $$ $$ P_3(x) = \prod_{i \in B_k} \frac{x-x_i}{x_3- x_i} = \frac{x - x_0}{x_3-x_0}\times \frac{x - x_1}{x_3-x_1}\times \frac{x - x_2}{x_3-x_2} $$ $f(x)$是一个过4个点的唯一的三次多项式。下面来验证一下拉格朗日多项式，假设有四个点(1,4)、(2,15)、(3,40)、(4,85)那么： $$ P_0(x) \times y_0= \frac{x - 2}{1 - 2}\times \frac{x - 3}{1 - 3} \times \frac{x - 4}{1- 4} \times 4= - \frac{(x -2)(x-3)(x-4)}{6} \times 4 = -\frac{4\times (x^3 -9x ^2 +26x -24)}{6} $$

$$ P_1(x) \times y_1= \frac{x - 1}{2 - 1}\times \frac{x - 3}{2 - 3} \times \frac{x - 4}{2- 4} \times 15= \frac{(x -1)(x-3)(x-4)}{2}\times 15 = \frac{45\times (x^3 - 8x^2 + 19x -12)}{6} $$

$$ P_2(x) \times y_2= \frac{x - 1}{3 - 1}\times \frac{x - 2}{3- 2} \times \frac{x - 4}{3- 4} \times 40= - \frac{(x -1)(x-2)(x-4)}{2}\times 40 = -\frac{120 \times(x^3 -7x^2 + 14x -8)}{6} $$

$$ P_3(x) \times y_3= \frac{x - 1}{4 - 1}\times \frac{x - 2}{4 - 2} \times \frac{x - 3}{4- 3} \times 85= \frac{(x -1)(x-2)(x-3)}{6}\times 85= \frac{85 \times(x^3 -6x^2 +11x-6)}{6} $$ $$\Longrightarrow f(x) = \frac{(-4 + 45 - 120 + 85)x^3}{6} + \frac{(36 -360 +840 -510)x^2} {6} + \frac{(-104 + 855 - 1680 + 935)x}{6} + \frac{96-540+960-510}{6} = x^3 + x^2 + x + 1$$

最后得到的$\;f(x) = x^3 + x^2 + x + 1$，可以代入$\;x = 4$，则有$\;f(4) = 64 + 16 + 4 + 1 = 85$。

11.2 SciPy拉格朗日插值

有了上边的理论不难理解拉格朗日插值多项式的意图了，编写一个简单的程序验证一下scipy的拉格朗日插值函数即scipy.interpolate模块下的 lagrange函数。

from scipy.interpolate  import lagrange
x = [1, 2, 3,4]
y = [4, 15, 40, 85]
ret = lagrange(x,y)
print ret

程序执行结果：

   3     2
1 x + 1 x + 1 x + 1

如果您当前的scipy不是1.10版，可以通过下面的命令升级scipy。 sudo pip install --upgrade scipy

感谢Klang(金浪)智能数据看板klang.org.cn鼎力支持！