- 拉格朗日插值
import matplotlib.pyplot as plt
from pylab import mpl
#计算插值多项式的系数
x = [0.2, 0.4, 0.6, 0.8, 1.0]
y = [0.98, 0.92, 0.81, 0.64, 0.38]
def ParametersOfLagrangeInterpolation(data_x,data_y,size):
parameters=[]
i=0;#i用来控制参数的个数
while i < size:
j = 0;#j用来控制循环的变量做累乘
temp = 1;
while j < size:
if(i != j):
temp*=data_x[i]-data_x[j]
j+=1;
parameters.append(data_y[i]/temp)
i += 1;
return parameters
#计算拉格朗日插值公式的值
def CalculateTheValueOfLarangeInterpolation(data_x,parameters,x):
returnValue=0
i = 0;
while i < len(parameters):
temp = 1
j = 0;
while j< len(parameters):
if(i!=j):
temp *=x-data_x[j]
j+=1
returnValue += temp * parameters[i]
i += 1
return returnValue
#将函数绘制成图像
def Draw(data_x,data_y,new_data_x,new_data_y):
plt.plot(new_data_x, new_data_y, label="拟合曲线", color="red")
plt.scatter(data_x,data_y, label="离散数据",color="yellow")
plt.scatter(1.75, 9.37890625, label="真实数据", color="orange")
plt.scatter(1.25, 2.44140625, color="green")
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
plt.title("Lagrange插值拟合数据")
plt.legend(loc="upper left")
plt.show()
parameters=ParametersOfLagrangeInterpolation(x,y,5)
datax=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2]
datay=[]
for temp in datax:
datay.append(CalculateTheValueOfLarangeInterpolation(x,parameters,temp))
x.append(1.75)
y.append(CalculateTheValueOfLarangeInterpolation(x,parameters,1.75))
Draw(x,y,datax,datay)
print("得到的四次Lagrange插值多项式为：L(x) = %f(x-0)(x-1)(x-1.5)(x-2) + %f(x-0)(x-0.5)(x-1.5)(x-2) + %f(x-0)(x-0.5)(x-1)(x-2) + %f(x-0)(x-0.5)(x-1)(x-1.5)"%(parameters[1],parameters[2],parameters[3],parameters[4]))
2.牛顿插值
import matplotlib.pyplot as plt
from pylab import mpl
import math
x = [0.2, 0.4, 0.6, 0.8, 1.0]
y = [0.98, 0.92, 0.81, 0.64, 0.38]
"""计算四次差商的值"""
def four_order_difference_quotient(x, y):
# i记录计算差商的次数，这里循4次，计算4次差商。
i = 0
quotient = [0, 0, 0, 0, 0]
while i < 4:
j = 4
while j > i:
if i == 0:
quotient[j]=((y[j]-y[j-1])/(x[j]-x[j-1]))
else:
quotient[j] = (quotient[j]-quotient[j-1])/(x[j]-x[j-1-i])
j -= 1
i += 1
return quotient;
def function(data):
return parameters[0] + parameters[1]*(data-0) + parameters[2]*(data-0)*(data-0.5) + parameters[3]*(data-0)*(data-0.5)*(data-1) + parameters[4]*(data-0)*(data-0.5)*(data-1)*(data-1.5)
"""计算插值多项式的值"""
def calculate_data(x,parameters):
returnData=[];
for data in x:
returnData.append(function(data))
return returnData
"""画函数的图像
newData为曲线拟合后的曲线
"""
def draw(newData):
plt.scatter(x,y,label="离散数据",color="blue")
plt.plot(datax,newData,label="牛顿插值拟合数据",color="orange")
plt.scatter(1.75,1.75**4, label="准确值",color="red")
plt.title("牛顿插值法拟合数据")
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
plt.legend(loc="upper left")
plt.show()
parameters=four_order_difference_quotient(x,y)
datax=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2]
datay=calculate_data(datax,parameters)
draw(datay)
print("牛顿插值多项式为：N(x) = %f(x-0) + %f(x-0)(x-0.5) + %f(x-0)(x-0.5)(x-1) + %f(x-0)(x-0.5)(x-1)(x-1.5)"%(parameters[1], parameters[2], parameters[3],parameters[4]))
3.三次样条插值
import numpy as np
import matplotlib.pyplot as plt
from pylab import mpl
"""
三次样条实现：
函数的自变量x:0.2, 0.4, 0.6, 0.8, 1.0
函数的因变量y:0.98, 0.92, 0.81, 0.64, 0.38
"""
x = [0.2, 0.4, 0.6, 0.8, 1.0]
y = [0.98, 0.92, 0.81, 0.64, 0.38]
"""
功能：完后对三次样条函数求解方程参数的输入
参数：要进行三次样条曲线计算的自变量
返回值：方程的参数
"""
def calculateEquationParameters(x):
'''
代解未知数为a1,b1,c1,d1,a2,b2,c2,d2,a3,b3,c3,d3
:param x:
:return: parameter
'''
# parameter为二维数组，用来存放参数，sizeOfInterval是用来存放区间的个数
parameter = []
n = len(x)
sizeOfInterval = n - 1 # 线段数
i = 1
# 首先输入方程两边相邻 节点 处函数值相等的方程为2(n-2)=4个方程; n个点要减去两个端点，然后乘以2
while i <= n-2:
# x[i]代入点i前一条线段的参数
data = np.zeros(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
# x[i]代入点i后一条线段的参数
data1 = np.zeros(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
parameter.append(data)
parameter.append(data1)
i += 1
# 输入 端点 处的函数值。为2个方程
data = np.zeros(sizeOfInterval * 4)
# 第一条线段
data[0] = x[0] * x[0] * x[0]
data[1] = x[0] * x[0]
data[2] = x[0]
data[3] = 1
parameter.append(data)
data = np.zeros(sizeOfInterval * 4)
# 最后一条线段
data[-4] = x[-1] * x[-1] * x[-1]
data[-3] = x[-1] * x[-1]
data[-2] = x[-1]
data[-1] = 1
parameter.append(data)
# 节点点函数一阶导数值相等为n-2=2个方程。
i = 1
while i <= n-2:
data = np.zeros(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)
parameter.append(data)
i += 1
# 节点函数二阶导数值相等为n-2=2个方程。且端点处的函数值的二阶导数为零，为2个方程。
i = 1
while i <= n-2:
data = np.zeros(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)
parameter.append(data)
i += 1
# 总共2(n-1)-2=10个方程
parameter = np.array(parameter)
return parameter[:, 2:] # 去掉前两个a1,b1
"""
功能：计算样条函数的系数。
参数：parametes为方程的系数，y为要插值函数的因变量。
返回值：三次插值函数的系数。
"""
def solutionOfEquation(parametes, y):
n = len(x)
sizeOfInterval = n - 1
result = np.zeros(sizeOfInterval * 4 - 2)
i = 1
# 节点处方程右边
while i < sizeOfInterval: # result[0,1,2,3]
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(parametes)
b = np.array(result)
return np.linalg.solve(a, b) # 解线性方程组
"""
功能：根据所给参数，计算三次函数的函数值：
参数:parameters为二次函数的系数，x为自变量
返回值：为函数的因变量
"""
def calculate(paremeters, x):
result = []
for data_x in x:
y = paremeters[0] * data_x * data_x * data_x + paremeters[1] * data_x * data_x + paremeters[2] * data_x + paremeters[3]
result.append(y)
return result
"""
功能：采点
参数：x
返回值：采取点和采取点的函数值
"""
def grasp_sample(x):
n = len(x)
i = 1
# result 为求解出来后的a1,b1,c1,d1,a2,b2,c2,d2,a3,b3,c3,d3
result = [0, 0]
temp = solutionOfEquation(calculateEquationParameters(x), y)
result.extend(temp)
samples_x = []
samples_y = []
# n-1段曲线
while i < n:
sample_x = np.arange(x[i-1], x[i], 0.01)
sample_y = calculate(result[(i-1)*4:i*4], sample_x)
samples_x.extend(sample_x)
samples_y.extend(sample_y)
i = i+1
samples_x.append(x[n-1])
samples_y.extend(calculate(result[-4:], [x[n-1]]))
return [samples_x, samples_y]
"""
功能：将函数绘制成图像
参数：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()
samples = grasp_sample(x)
# print(samples[0])
# print(samples[1])
Draw(x, y, samples[0], samples[1])
3.
import matplotlib.pyplot as plt
from scipy import interpolate
import numpy as np
import matplotlib.font_manager as mpt
x =[0,1,4,9,16,25,36,49,64]#创建x列表储存x数据
y =[0,1,2,3,4,5,6,7,8]#创建y列表储存y数据
# 进行样条插值
tck = interpolate.splrep(x,y)
xx = np.linspace(min(x),max(x),100)
yy = interpolate.splev(xx,tck,der=0)
print(xx)
# 画图
plt.plot(x,y,'o',xx,yy)
plt.legend(['true','Cubic-Spline'])
plt.xlabel('x', fontproperties="zhfont") #注意后面的字体属性
plt.ylabel('y')
plt.title('S(X)', fontproperties="zhfont")
# 保存图片
plt.savefig('out.jpg')
plt.show()
扫一扫
4万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
