实现拉格朗日插值方法 (Python)

chloeee129
·
(修改过)
·
IPFS
我学python还没到人生苦短的阶段,还是嫌命长

 【还没研究过blogspot的排版方式,好像不能插入latex?所以先用截图来展示表达式了】

有高斯形状的分布函数:


其中\mu = 50, \sigma = 15令x0 = 5.0, xn = 95.0, 在[x0, xn]范围内,等间距地选取 n+1 个值

以下为n = 2, 4, 6 时的多项式插值函数和原高斯函数的绘制代码:

#对高斯函数做拉格朗日插值近似

import numpy as np

import matplotlib.pyplot as plt

plt.rcParams['xtick.labelsize'] = 20

plt.rcParams['ytick.labelsize'] = 20

#绘制n=2,4,6时的多项式差值函数

x0 = 5.0

xn = 95.0

for n in [2, 4, 6]:

   x_arr = np.arange(x0, xn+(xn-x0)/n, (xn-x0)/n)

   y_arr = np.exp(-(x_arr-50)**2/(2*15**2))/(15*np.sqrt(2*np.pi))

   x = np.arange(x0, xn+0.01, 0.01)

   lk = [y_arr[i] for i in range(len(x_arr))]

   for i in range(len(x_arr)):

       for j in range(n+1):

           if j == i:

               continue

           else:

               lk[i] *= (x-x_arr[j])/(x_arr[i]-x_arr[j])

   Pn = sum(lk)

   plt.figure(1,figsize=(15,8),dpi=300)

   lbl='n='+str(n)

   plt.plot(x,Pn,label=lbl)

   plt.plot(x_arr,y_arr,'ko',ms=10)


#绘制高斯函数

x_arr = np.arange(x0, xn+0.1, 0.1)

y_arr = np.exp(-(x_arr-50)**2/(2*15**2))/(15*np.sqrt(2*np.pi))

plt.plot(x_arr,y_arr,label='Guass')

plt.legend(prop={'size': 16})

plt.grid()

______________________________

下图为output:


 

 

 下面附赠线性插值的画法,可以参考如何用python写分段函数。

代码:

import numpy as np

import matplotlib.pyplot as plt

plt.rcParams['xtick.labelsize'] = 20

plt.rcParams['ytick.labelsize'] = 20


#用线性插值法近似高斯函数

x0 = 5.0

xn = 95.0

for n in [2, 4, 6]:

   x_arr = np.arange(x0, xn+(xn-x0)/n, (xn-x0)/n)

   y_arr = np.exp(-(x_arr-50)**2/(2*15**2))/(15*np.sqrt(2*np.pi))

   #定义分段函数

   def fx(x):

       for j in range(n):

           if x_arr[j] <= x <= x_arr[j+1]:

               return y_arr[j]*(x-x_arr[j+1])/(x_arr[j]-x_arr[j+1])+y_arr[j+1]*(x-x_arr[j])/(x_arr[j+1]-x_arr[j])

           else:

               continue

   x = np.arange(x0, xn+1, 1)

   y = np.array([fx(t) for t in x])

   plt.figure(1,figsize=(15,8),dpi=300)

   lbl='n='+str(n)

   plt.plot(x,y,label=lbl)

   plt.plot(x_arr,y_arr,'ko',ms=10)

plt.legend(prop={'size': 16})

plt.grid()

______________________________

 

 下图为输出:



 

 

 以上禁止转载,仅供参考@奇特比喻

CC BY-NC-ND 2.0 授权

喜欢我的作品吗?别忘了给予支持与赞赏,让我知道在创作的路上有你陪伴,一起延续这份热忱!