拉格朗日插值与的Python实现

最新推荐文章于 2024-07-06 16:48:32 发布

原创最新推荐文章于 2024-07-06 16:48:32 发布 · 679 阅读

2 ·

CC 4.0 BY-SA版权

问题解决专栏收录该内容

5 篇文章

订阅专栏

本文深入探讨了拉格朗日插值法的基本原理，包括线性插值和高次插值的数学推导过程，展示了如何利用Python实现拉格朗日插值，并通过实例比较了不同插值方法的效果。

拉格朗日插值数学原理

线性插值

对于一个函数 $f (x)$ 在区间 $[a, b]$ 上有定义，且已知在点 $a<=x_0<x_1<···<x_n<=b$ 上的值 $y_0,y_1,···,y_n$ . 若存在一个函数 $P(x)=a_0+a_1x+···+a_nx$ .使得：
$P(x_i) = y_i$
我们就说 $P (x) 为 f (x)$ 的插值函数。对于如何求出这个插值函数，我们有多种方法，我们先设 $n = 1$ 时的简单情况，假定给定区间 $x_k,x_{k+1}]$ ,及插值端点值 $y_k=f(x_k),y_{k+1}=f(x_{k+1})$ ,要求线性多项式 $L_1(x)$ 满足：
$L_1(x_k)=y_k,L_1(x_{k+1})=y_{k+1}$
上式的几何意义是通过两点 $x_k,y_k),(x_{k+1},y_{k+1})$ 的直线，用公式表达它们的几何意义：
$L_1(x) = y_k+\frac{y_{k+1}-y_k}{x_{k+1}-x_k}(x-x_k) (点斜式) \\ L_1(x) = y_k\frac{x_{k+1}-x}{x_{k_1}-x_k}+y_{k+1}\frac{x-x_k}{x_{k+1}-x_k} (两点式)$
由以上两式可以看出 $L_1(1x)$ 是由两个线性函数：
$l_k(x)=\frac{y_{k+1}-y_k}{x_{k+1}-x_k}，l_{k+1}=\frac{x-x_k}{x_{k+1}-x_k}$
线性组合得到的，其系数分别为 $y_k$ 及 $y_{k+1}$ ,即：
$L_1(x)=y_1l_k(x)+y_{k+1}l_{k+1}(x).$
我们称 $l_k(x)$ 及 $k_{k+1}$ 为线性插值基函数， $L (x)$ 为一次插值多项式，同理当 $n = 2$ 时我们也可以构造出二次插值多项式，此外我们可以把它推广到 $n$ ，所以就得到了 $n$ 次插值多项式，如下公式；
$l_k(x) = \frac{(x-x_0)···(x-x_{k-1})(x-x_{k+1})···(x-x_{n})}{(x_k-x_0)···(x_k-x_{k-1})(x_k-x_{k+1})···(x_k-x_{n})} \\ L_n(x)=\sum_{k=1}^ny_kl_k(x)\\ L_n(x_j)\sum_{k=1}^ny_kl_k(x_j) y_j , j= 0,1,···,n. (1)$
我们得到的(1)式被称为拉格朗日插值多项式。

算法实现

import matplotlib.pyplot as plt
import numpy as np
from numpy import poly1d
from scipy._lib.six import xrange
# 自己实现的拉格朗日插值
def h(x,y,a):
    if len(x)!= len(y):
        print('向量长度必须相等！')
        exit(-1)
    ans=0.0
    for i in range(len(y)):
        temp = 1.0
        t = y[i]
        #这个循环算出一个插值基函数值
        for j in range(len(y)):
            if i !=j:
                temp *= (a-x[j])/(x[i]-x[j])
        t *= temp #得出一个插值多项式
        ans +=t
    return ans

def lagrange(x,w):
    M = len(x)
    p = poly1d(0.0)
    for j in xrange(M):
        pt = poly1d(w[j])
        print(pt)
        for k in xrange(M):
            if k == j:
                continue
            fac = x[j] - x[k]
            pt *= poly1d([1.0, -x[k]]) / fac
        p += pt
    return p

origin_x = np.linspace(-4, 4, 50)
origin_y = fun(origin_x)

x = np.linspace(-4, 4, 9)
y = fun(x)
f = lagrange(x,y)

Inter_x = np.linspace(-4, 4, 50)

NormlagInter_y = [ f(a) for a in Inter_x ]
mylagInter_y = [neton_inter(x, y, a) for a in Inter_x]

plt.figure(figsize=[15, 10])
ax = plt.subplot(1, 3, 1)


ax.set_title("orgin")
plt.plot(origin_x, origin_y, c='red')

ax = plt.subplot(1,3,2)
ax.set_title("NormlagInter")
plt.plot(Inter_x,NormlagInter_y)

ax = plt.subplot(1,3,3)
ax.set_title('mylagInter')
plt.plot(Inter_x,mylagInter_y,c='blue')
plt.show()