拉格朗日插值与的Python实现

mac2024-05-22  34

拉格朗日插值数学原理

线性插值

对于一个函数 f ( x ) f(x) f(x)在区间 [ a , b ] [a,b] [a,b]上有定义,且已知在点 a < = x 0 < x 1 < ⋅ ⋅ ⋅ < x n < = b a<=x_0<x_1<···<x_n<=b a<=x0<x1<<xn<=b上的值 y 0 , y 1 , ⋅ ⋅ ⋅ , y n y_0,y_1,···,y_n y0,y1,,yn. 若存在一个函数 P ( x ) = a 0 + a 1 x + ⋅ ⋅ ⋅ + a n x P(x)=a_0+a_1x+···+a_nx P(x)=a0+a1x++anx.使得: P ( x i ) = y i P(x_i) = y_i P(xi)=yi 我们就说 P ( x ) 为 f ( x ) P(x)为f(x) P(x)f(x)的插值函数。对于如何求出这个插值函数,我们有多种方法,我们先设 n = 1 n=1 n=1时的简单情况,假定给定区间 [ x k , x k + 1 ] [x_k,x_{k+1}] [xk,xk+1],及插值端点值 y k = f ( x k ) , y k + 1 = f ( x k + 1 ) y_k=f(x_k),y_{k+1}=f(x_{k+1}) yk=f(xk),yk+1=f(xk+1),要求线性多项式 L 1 ( x ) L_1(x) L1(x)满足: L 1 ( x k ) = y k , L 1 ( x k + 1 ) = y k + 1 L_1(x_k)=y_k,L_1(x_{k+1})=y_{k+1} L1(xk)=yk,L1(xk+1)=yk+1 上式的几何意义是通过两点 ( x k , y k ) , ( x k + 1 , y k + 1 ) (x_k,y_k),(x_{k+1},y_{k+1}) (xk,yk),(xk+1,yk+1)的直线,用公式表达它们的几何意义: L 1 ( x ) = y k + y k + 1 − y k x k + 1 − x k ( x − x k ) ( 点 斜 式 ) L 1 ( x ) = y k x k + 1 − x x k 1 − x k + y k + 1 x − x k x k + 1 − x k ( 两 点 式 ) 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} (两点式) L1(x)=yk+xk+1xkyk+1yk(xxk)()L1(x)=ykxk1xkxk+1x+yk+1xk+1xkxxk() 由以上两式可以看出 L 1 ( 1 x ) L_1(1x) L1(1x)是由两个线性函数: l k ( x ) = y k + 1 − y k x k + 1 − x k , l k + 1 = x − x k x k + 1 − x k 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} lk(x)=xk+1xkyk+1yklk+1=xk+1xkxxk 线性组合得到的,其系数分别为 y k y_k yk y k + 1 y_{k+1} yk+1,即: L 1 ( x ) = y 1 l k ( x ) + y k + 1 l k + 1 ( x ) . L_1(x)=y_1l_k(x)+y_{k+1}l_{k+1}(x). L1(x)=y1lk(x)+yk+1lk+1(x). 我们称 l k ( x ) l_k(x) lk(x) k k + 1 k_{k+1} kk+1为线性插值基函数, L ( x ) L(x) L(x)为一次插值多项式,同理当 n = 2 n=2 n=2时我们也可以构造出二次插值多项式,此外我们可以把它推广到 n n n,所以就得到了 n n n次插值多项式,如下公式; l k ( x ) = ( 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 ) = ∑ k = 1 n y k l k ( x ) L n ( x j ) ∑ k = 1 n y k l k ( x j ) y j , j = 0 , 1 , ⋅ ⋅ ⋅ , n . ( 1 ) 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) lk(x)=(xkx0)(xkxk1)(xkxk+1)(xkxn)(xx0)(xxk1)(xxk+1)(xxn)Ln(x)=k=1nyklk(x)Ln(xj)k=1nyklk(xj)yj,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()
最新回复(0)