python分段线性插值_计算⽅法(3)——分段插值法(附
Python程序)
在上⼀节计算⽅法(2)——插值法(附Python程序)当中,主要讲了插值法,介绍了龙格现象,并给出了插值法的代码。
这⼀讲主要分段插值中的分段线性插值和分段Hermite插值,并给出分段插值的Python程序。
在此之前需要注意⼀下,n为区间数,n+1为插值节点的个数。
分段线性插值
分段线性插值,需要两个列表,⼀个⽤于存放各点的x坐标,⼀个⽤于存放各点的y坐标。
因为分段插值的算法需要x坐标按顺序增长,⽽调⽤该函数时可能传⼊的是⼀些未经排序的散点,所以在传⼊xlist和ylist时,需要对xlist进⾏排序,同时ylist也要跟着变化。解决⽅法是转换成字典,再对字典的key值排序。
# 分段线性插值,需要两个列表
def segmented_linear_interpolate(xlist,ylist,x):
"""n = # of intervals, which is derived from len of xlistlen of xlist is always one item biger than # of intervals"""
# we have to make sure that items in xlist is in order
data = dict(zip(xlist,ylist))
# 按照key排序,也就是xlist
data = sorted(data.items(),key=lambda item:item[0])
data = dict(data)
xlist = list(data.keys())
ylist = list(data.values())
n = len(xlist)-1
if n == 0:
raise ValueError("n should be greater or equal to 1")
# print("segmented interpolate, n =",n)
# 需要把新来的元素判断⼀下在哪个区间
i = -1
for t in xlist:
if x >= t:
i += 1
if i == -1 or i > len(xlist)-1:
raise ValueError("x should be between%fand%f"%(xlist[0],xlist[-1]))
if i == len(xlist)-1:
return ylist[i]
return (x-xlist[i+1])/(xlist[i]-xlist[i+1])*ylist[i] + (x-xlist[i])/(xlist[i+1]-xlist[i])*ylist[i+1]
程序⽐较臃肿,还有优化空间。下⾯看⼀下效果:
if __name__ == "__main__":
import numpy as np
import matplotlib.pyplot as plt
f = lambda x: 1/(1+25*x**2)
# 待插值的元素值
x_points = [0,1,2,3,4]
y_points = [20,18,10,2,1]
# 分段线性插值
f = lambda t: segmented_linear_interpolate(x_points,y_points,t)
x = np.linspace(0,4)
y = list(map(f,x))
# 画图
plt.figure("segmented linear interpolation")
plt.scatter(x_points,y_points,color = "orange")
plt.plot(x,y)
plt.legend(["segmented linear interpolation","scattered points"])
plt.show()
这⼀段程序,根据传⼊的五个数据点,⾃动分了四段进⾏插值,所得函数如图。分段线性插值(n=4)
分段线性插值可以很好地解决上⼀节所讲的龙格现象。上⼀节讲的龙格现象(拉格朗⽇插值 n=10)龙格现象的解决(分段线性插值 n=10)⽤分段线性插值解决龙格现象的代码:
# 龙格现象的解决-分段线性插值
if __name__ == "__main__":
import numpy as np
import matplotlib.pyplot as plt
f = lambda x: 1/(1+25*x**2)
# 待插值的元素值
x_points = np.linspace(-1,1,11)
y_points = list(map(f,x_points))
# 分段线性插值
fx = lambda t: segmented_linear_interpolate(x_points,y_points,t)
x = np.linspace(-1,1,51)
y = list(map(fx,x))
# 画图
plt.figure("segmented interpolation")
plt.scatter(x_points,y_points,color = "orange")
plt.plot(x,y)
plt.legend(["segmented interpolation","scattered points"])
plt.show()
上⼀节⽤拉格朗⽇插值产⽣龙格现象的代码:
# 龙格现象的产⽣
if __name__ == "__main__":
import numpy as np
import matplotlib.pyplot as plt
f = lambda x: 1/(1+25*x**2)
# 待插值的元素值
x_points = np.linspace(-1,1,11)
print(x_points)
y_points = list(map(f,x_points))
# 拉格朗⽇插值
x = np.linspace(-1,1)
y = list(map(lambda t: lagrange_interpolate(x_points,y_points,t),x))
# 画图
plt.figure("lagrange interpolation")
plt.scatter(x_points,y_points,color = "orange")
plt.plot(x,y)
plt.legend(["lagrange interpolation","scattered points"])
plt.show()
分段Hermite插值
分段Hermite插值,需要三个列表存放数据:各点的x坐标,各点的y坐标,各点的导数值。
# 分段Hermite插值,需要输⼊三个列表
linspace numpydef segmented_hermite_interpolate(x_list,y_list,y_prime_list,x):
"""n = # of intervals, which is derived from len of xlistlen of xlist is always one item biger than # of intervals""" #we have to make sure that items in xlist is in order
# 按照x_list给y_list排序
data = dict(zip(x_list,y_list))
data = sorted(data.items(),key=lambda item:item[0])
data = dict(data)
xlist = list(data.keys())
ylist = list(data.values())
# 按照x_list给y_prime_list排序
data = dict(zip(x_list,y_prime_list))
data = sorted(data.items(),key=lambda item:item[0])
data = dict(data)
y_prime_list = list(data.values())
n = len(xlist)-1
if n == 0:
raise ValueError("n should be greater or equal to 1")
# print("segmented interpolate, n =",n)
# 需要把新来的元素判断⼀下在哪个区间
i = -1
for t in xlist:
if x >= t:
i += 1
if i == -1 or i > len(xlist)-1:
raise ValueError("x should be between%fand%f"%(xlist[0],xlist[-1]))
if i == len(xlist)-1:
return ylist[i]
alpha0 = lambda x: ((x-xlist[i+1])/(xlist[i]-xlist[i+1]))**2 * (2*(x-xlist[i])/(xlist[i+1]-xlist[i])+1) alpha1 = lambda x: ((x-xlist[i])/(xlist[i+1]-xlist[i]))**2 * (2*(x-xlist[i+1])/(xlist[i]-xlist[i+1])+1)
beta0 = lambda x: ((x-xlist[i+1])/(xlist[i]-xlist[i+1]))**2 * (x-xlist[i])
beta1 = lambda x: ((x-xlist[i])/(xlist[i+1]-xlist[i]))**2 * (x-xlist[i+1])
H = alpha0(x)*ylist[i] + alpha1(x)*ylist[i+1] + beta0(x)*y_prime_list[i] + beta1(x)*y_prime_list[i+1] return H
同样还有优化的空间,⽐如排序的时候可以同时根据x_list来排y_list和y_prime_list,⽽不是排两次。
运⾏⼀下程序:
if __name__ == "__main__":
import numpy as np
import matplotlib.pyplot as plt
# 待插值的元素值
x_points = [0,1,2,3,4]
y_points = [5,4,3,2,1]
y_primes = [0,0,0,0,0] # 让各点导数等于0
# 分段hermite插值
f = lambda t: segmented_hermite_interpolate(x_points,y_points,y_primes,t)
x = np.linspace(0,4)
y = list(map(f,x))
# 画图
plt.figure("segmented hermite interpolation")
plt.scatter(x_points,y_points,color = "orange")
plt.plot(x,y)
plt.legend(["segmented hermite interpolation","scattered points"])
plt.show()分段Hermite插值的效果
可以看到分段Hermite插值,不仅满⾜了各点的位置要求,还满⾜了各点导数值等于0的要求。-END-

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。