https://blog.csdn.net/qq_36607894/article/details/100012334?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522162678321416780269872462%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=162678321416780269872462&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~top_positive~default-1-100012334.first_rank_v2_pc_rank_v29&utm_term=%E6%8F%92%E5%80%BC&spm=1018.2226.3001.4187
一、多项式插值
二、分段线性插值
#Program 0.6 Linear Interploation
import numpy as np
import matplotlib.pyplot as plt
#分段线性插值闭包
def get_line(xn, yn):
def line(x):
index = -1
#找出x所在的区间
for i in range(1, len(xn)):
if x <= xn[i]:
index = i-1
break
else:
i += 1
if index == -1:
return -100
#插值
result = (x-xn[index+1])*yn[index]/float((xn[index]-xn[index+1])) + (x-xn[index])*yn[index+1]/float((xn[index+1]-xn[index]))
return result
return line
xn = [i for i in range(-50,50,10)]
yn = [i**2 for i in xn]
#分段线性插值函数
lin = get_line(xn, yn)
x = [i for i in range(-50, 40)]
y = [lin(i) for i in x]
#画图
plt.plot(xn, yn, 'ro')
plt.plot(x, y, 'b-')
plt.show()三、Hermite插值
不但要求在节点上的函数值相等，而且还要求对应的导数值也相等,甚至要求高阶导数也相等，满足这种要求的插值多项式就是埃尔米特插值多项式。
# -*- coding: utf-8 -*-
#Program 0.5 Hermite Interpolation
import matplotlib.pyplot as plt
import numpy as np
#计算基函数的导数值
def dl(i, xi):
result = 0.0
for j in range(0,len(xi)):
if j!=i:
result += 1/(xi[i]-xi[j])
result *= 2
return result
#计算基函数值
def l(i, xi, x):
deno = 1.0
nu = 1.0
for j in range(0, len(xi)):
if j!= i:
deno *= (xi[i]-xi[j])
nu *= (x-xi[j])
return nu/deno
#Hermite插值函数
def get_Hermite(xi, yi, dyi):
def he(x):
result = 0.0
for i in range(0, len(xi)):
result += (yi[i]+(x-xi[i])*(dyi[i]-2*yi[i]*dl(i, xi))) * ((l(i,xi,x))**2)
return result
return he
import math
sr_x = [(i * math.pi) + (math.pi / 2) for i in range(-3, 3)]
sr_fx = [math.sin(i) for i in sr_x]
deriv = [0 for i in sr_x] # 导数都为 0
Hx = get_Hermite(sr_x, sr_fx, deriv) # 获得插值函数
tmp_x = [i * 0.1 * math.pi for i in range(-20, 20)] # 测试用例
tmp_y = [Hx(i) for i in tmp_x] # 根据插值函数获得测试用例的纵坐标
#画图
plt.plot(sr_x, sr_fx, 'ro')
plt.plot(tmp_x, tmp_y, 'b-')
plt.title('Hermite Interpolation')
plt.show()四、样条插值
1.2.2 三次样条插值
定义
每个子区间的多项式函数是三次的，所以名字里有“三次”，就像前面拉格朗日插值的分段二次插值的每个子区间的插值函数是二次的一样
整个插值函数曲线由分段三次曲线并接而成，并且在连接点也就是样点上必须要二阶连续可导
n+1个数据点，共有n个分段多项式
相邻两个子多项式之间必须在函数值，一阶导数，二阶导数上都相等
对比发现，分段三次hermite插值只考虑了一阶导数，分段三次样条考虑了一阶导数和二阶导数，所以后者得到的结果更加精确
import numpy as np
import matplotlib.pyplot as plt
from pylab import mpl
"""
三次样条实现：
函数x的自变量为:3, 4.5, 7, 9
因变量为：2.5, 1 2.5, 0.5
"""
x = [3, 4.5, 7, 9]
y = [2.5, 1, 2.5, 0.5]
"""
功能：完后对三次样条函数求解方程参数的输入
参数：要进行三次样条曲线计算的自变量
返回值：方程的参数
"""
def calculateEquationParameters(x):
#parameter为二维数组，用来存放参数，sizeOfInterval是用来存放区间的个数
parameter = []
sizeOfInterval=len(x)-1;
i = 1
#首先输入方程两边相邻节点处函数值相等的方程为2n-2个方程
while i < len(x)-1:
data = init(sizeOfInterval*4)
data[(i-1)*4] = x[i]*x[i]*x[i]
data[(i-1)*4+1] = x[i]*x[i]
data[(i-1)*4+2] = x[i]
data[(i-1)*4+3] = 1
data1 =init(sizeOfInterval*4)
data1[i*4] =x[i]*x[i]*x[i]
data1[i*4+1] =x[i]*x[i]
data1[i*4+2] =x[i]
data1[i*4+3] = 1
temp = data[2:]
parameter.append(temp)
temp = data1[2:]
parameter.append(temp)
i += 1
# 输入端点处的函数值。为两个方程, 加上前面的2n - 2个方程，一共2n个方程
data = init(sizeOfInterval * 4 - 2)
data[0] = x[0]
data[1] = 1
parameter.append(data)
data = init(sizeOfInterval * 4)
data[(sizeOfInterval - 1) * 4 ] = x[-1] * x[-1] * x[-1]
data[(sizeOfInterval - 1) * 4 + 1] = x[-1] * x[-1]
data[(sizeOfInterval - 1) * 4 + 2] = x[-1]
data[(sizeOfInterval - 1) * 4 + 3] = 1
temp = data[2:]
parameter.append(temp)
# 端点函数一阶导数值相等为n-1个方程。加上前面的方程为3n-1个方程。
i=1
while i < sizeOfInterval:
data = init(sizeOfInterval * 4)
data[(i - 1) * 4] = 3 * x[i] * x[i]
data[(i - 1) * 4 + 1] = 2 * x[i]
data[(i - 1) * 4 + 2] = 1
data[i * 4] = -3 * x[i] * x[i]
data[i * 4 + 1] = -2 * x[i]
data[i * 4 + 2] = -1
temp = data[2:]
parameter.append(temp)
i += 1
# 端点函数二阶导数值相等为n-1个方程。加上前面的方程为4n-2个方程。且端点处的函数值的二阶导数为零，为两个方程。总共为4n个方程。
i = 1
while i < len(x) - 1:
data = init(sizeOfInterval * 4)
data[(i - 1) * 4] = 6 * x[i]
data[(i - 1) * 4 + 1] = 2
data[i * 4] = -6 * x[i]
data[i * 4 + 1] = -2
temp = data[2:]
parameter.append(temp)
i += 1
return parameter
"""
对一个size大小的元组初始化为0
"""
def init(size):
j = 0;
data = []
while j < size:
data.append(0)
j += 1
return data
"""
功能：计算样条函数的系数。
参数：parametes为方程的系数，y为要插值函数的因变量。
返回值：三次插值函数的系数。
"""
def solutionOfEquation(parametes,y):
sizeOfInterval = len(x) - 1;
result = init(sizeOfInterval*4-2)
i=1
while i<sizeOfInterval:
result[(i-1)*2]=y[i]
result[(i-1)*2+1]=y[i]
i+=1
result[(sizeOfInterval-1)*2]=y[0]
result[(sizeOfInterval-1)*2+1]=y[-1]
a = np.array(calculateEquationParameters(x))
b = np.array(result)
for data_x in b:
print(data_x)
return np.linalg.solve(a,b)
"""
功能：根据所给参数，计算三次函数的函数值：
参数:parameters为二次函数的系数，x为自变量
返回值：为函数的因变量
"""
def calculate(paremeters,x):
result=[]
for data_x in x:
result.append(paremeters[0]*data_x*data_x*data_x+paremeters[1]*data_x*data_x+paremeters[2]*data_x+paremeters[3])
return result
"""
功能：将函数绘制成图像
参数：data_x,data_y为离散的点.new_data_x,new_data_y为由拉格朗日插值函数计算的值。x为函数的预测值。
返回值：空
"""
def Draw(data_x,data_y,new_data_x,new_data_y):
plt.plot(new_data_x, new_data_y, label="拟合曲线", color="black")
plt.scatter(data_x,data_y, label="离散数据",color="red")
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
plt.title("三次样条函数")
plt.legend(loc="upper left")
plt.show()
result=solutionOfEquation(calculateEquationParameters(x),y)
new_data_x1=np.arange(3, 4.5, 0.1)
new_data_y1=calculate([0,0,result[0],result[1]],new_data_x1)
new_data_x2=np.arange(4.5, 7, 0.1)
new_data_y2=calculate([result[2],result[3],result[4],result[5]],new_data_x2)
new_data_x3=np.arange(7, 9.5, 0.1)
new_data_y3=calculate([result[6],result[7],result[8],result[9]],new_data_x3)
new_data_x=[]
new_data_y=[]
new_data_x.extend(new_data_x1)
new_data_x.extend(new_data_x2)
new_data_x.extend(new_data_x3)
new_data_y.extend(new_data_y1)
new_data_y.extend(new_data_y2)
new_data_y.extend(new_data_y3)
Draw(x,y,new_data_x,new_data_y)
五、N维插值
一维插值：
# -*-coding:utf-8 -*-
import numpy as np
from scipy import interpolate
import pylab as pl
x=np.linspace(0,10,11)
#x=[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]
y=np.sin(x)
xnew=np.linspace(0,10,101)
pl.plot(x,y,"ro")
for kind in ["nearest","zero","slinear","quadratic","cubic"]:#插值方式
#"nearest","zero"为阶梯插值
#slinear 线性插值
#"quadratic","cubic" 为2阶、3阶B样条曲线插值
f=interpolate.interp1d(x,y,kind=kind)
# ‘slinear’, ‘quadratic’ and ‘cubic’ refer to a spline interpolation of first, second or third order)
ynew=f(xnew)
pl.plot(xnew,ynew,label=str(kind))
pl.legend(loc="lower right")
pl.show()
二维插值：
# -*- coding: utf-8 -*-
"""
演示二维插值。
"""
import numpy as np
from scipy import interpolate
import pylab as pl
import matplotlib as mpl
def func(x, y):
return (x+y)*np.exp(-5.0*(x**2 + y**2))
# X-Y轴分为15*15的网格
y,x= np.mgrid[-1:1:15j, -1:1:15j]
fvals = func(x,y) # 计算每个网格点上的函数值 15*15的值
print len(fvals[0])
#三次样条二维插值
newfunc = interpolate.interp2d(x, y, fvals, kind='cubic')
# 计算100*100的网格上的插值
xnew = np.linspace(-1,1,100)#x
ynew = np.linspace(-1,1,100)#y
fnew = newfunc(xnew, ynew)#仅仅是y值 100*100的值
# 绘图
# 为了更明显地比较插值前后的区别，使用关键字参数interpolation='nearest'
# 关闭imshow()内置的插值运算。
pl.subplot(121)
im1=pl.imshow(fvals, extent=[-1,1,-1,1], cmap=mpl.cm.hot, interpolation='nearest', origin="lower")#pl.cm.jet
#extent=[-1,1,-1,1]为x,y范围 favals为
pl.colorbar(im1)
pl.subplot(122)
im2=pl.imshow(fnew, extent=[-1,1,-1,1], cmap=mpl.cm.hot, interpolation='nearest', origin="lower")
pl.colorbar(im2)
pl.show()