一、问题引入
对于给出如下的离散的数据点,现在想根据如下的数据点来推测x=5时的值,我们应该采用什么方法呢?
x | f ( x) |
3.0 | 2.5 |
4.5 | 1.0 |
7.0 | 2.5 |
9.0 | 0.5 |
我们知道在平面上两个点确定一条直线,三个点确定一条抛物线(假设曲线的类型是抛物线),那么现在有四个点,我们很自然的会想到,既然两个点确定一条直线,那么最简单的方法就是,两个点之间连一条线,两个点之间连一条线,最后得到的一种折线图如下:这样我们只要确定x=5时的直线,把自变量的值带进去,就显然会得到预测的函数值。但是,这种方法显然具有很明显的缺陷:曲线不够光滑,连接点处的斜率变化过大。可能会导致函数的一阶导数不连续。那么我们应该如何解决这个问题呢?
二次样条的原理
我们会想到既然直线不行,那么我们就用曲线来近似的代替和描述。最简单的曲线是二次函数,如果我们用二次函数:aX^2+bx+c来描述曲线,最后的结果可能会好一点,下图中一共有4个点,可以分成3个区间。每一个区间都需要一个二次函数来描述,一共需要9个未知数。下面的任务就是找出9个方程。
如下图所示:一共有x0,x1,x2,x3四个点,三个区间,每个区间上都有一个方程。
1>曲线方程在节点处的值必须相等,即函数在x1,x2两个点处的值必须符合两个方程,这里一共是4个方程:
a1*x1^2+b1*x1+c1=f(x1)
a2*x1^2+b2*x1+c2=f(x1)
a2*x2^2+b2*x2+c2=f(x2)
a3*x2^2+b3*x2+c3=f(x2)
2>第一个端点和最后一个端点必须过第一个和最后一个方程:这里一共是2个方程
3>节点处的一阶导数的值必须相等。这里为两个方程。
2*a1*x1+b1=2*a2*x1+b2
2*a2*x2+b2=2*a3*x2+b3
4>在这里假设第一个方程的二阶导数为0:这里为一个方程:
a1=0
上面是对应的9个方程,现在只要把九个方程联立求解,最后就可以实现预测x=5处节点的值。
下面是写成矩阵的形式,由于a1=0,所以未知数的个数少了一个。
下面是二次样条的python实现,最后将结果绘制在图上。
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]
"""一共有三个区间,用二次样条求解,需要有9个方程"""
"""
功能:完后对二次样条函数求解方程参数的输入
参数:要进行二次样条曲线计算的自变量
返回值:方程的参数
"""
def calculateEquationParameters(x):
#parameter为二维数组,用来存放参数,sizeOfInterval是用来存放区间的个数
parameter = []
sizeOfInterval=len(x)-1;
i = 1
#首先输入方程两边相邻节点处函数值相等的方程为2n-2个方程
while i < len(x)-1:
data = init(sizeOfInterval*3)
data[(i-1)*3]=x[i]*x[i]
data[(i-1)*3+1]=x[i]
data[(i-1)*3+2]=1
data1 =init(sizeOfInterval*3)
data1[i * 3] = x[i] * x[i]
data1[i * 3 + 1] = x[i]
data1[i * 3 + 2] = 1
temp=data[1:]
parameter.append(temp)
temp=data1[1:]
parameter.append(temp)
i += 1
#输入端点处的函数值。为两个方程,加上前面的2n-2个方程,一共2n个方程
data = init(sizeOfInterval*3-1)
data[0] = x[0]
data[1] = 1
parameter.append(data)
data = init(sizeOfInterval *3)
data[(sizeOfInterval-1)*3+0] = x[-1] * x[-1]
data[(sizeOfInterval-1)*3+1] = x[-1]
data[(sizeOfInterval-1)*3+2] = 1
temp=data[1:]
parameter.append(temp)
#端点函数值相等为n-1个方程。加上前面的方程为3n-1个方程,最后一个方程为a1=0总共为3n个方程
i=1
while i < len(x) - 1:
data = init(sizeOfInterval * 3)
data[(i - 1) * 3] =2*x[i]
data[(i - 1) * 3 + 1] =1
data[i*3]=-2*x[i]
data[i*3+1]=-1
temp=data[1:]
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*3-1)
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)
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+paremeters[1]*data_x+paremeters[2])
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,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]],new_data_x2)
new_data_x3=np.arange(7, 9.5, 0.1)
new_data_y3=calculate([result[5],result[6],result[7]],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)
二次样条函数运行之后的结果如下,从图像中,我们可以看出,二次样条在函数的连接处的曲线是光滑的。这时候,我们将x=5输入到函对应的函数端中,就可以预测相应的函数值。但是,这里还有一个问题,就是二次样条函数的前两个点是直线,而且函数的最后一个区间内看起来函数凸出很高。我们还想解决这些问题,这时候,我们想是否可以用三次样条函数来进行函数的模拟呢?
三次样条的原理:
三次样条的原理和二次样条的原理相同,我们用函数aX^3+bX^2+cX+d这个函数来进行操作,这里一共是4个点,分为3个区间,每个区间一个三次样条函数的话,一共是12个方程,只要我们找出这12个方程,这个问题就算解决了。
1>内部节点处的函数值应该相等,这里一共是4个方程。
2>函数的第一个端点和最后一个端点,应该分别在第一个方程和最后一个方程中。这里是2个方程。
3>两个函数在节点处的一阶导数应该相等。这里是两个方程。
4>两个函数在节点处的二阶导数应该相等,这里是两个方程。
5>端点处的二阶导数为零,这里是两个方程。
a1=0
b1=0
端点条件
上面我们对端点处的样条进行了假设,为什么呢?其实端点可以有多种不同的限制,常见有3种。
自由边界 Natural
首尾两端没有受到任何使他们弯曲的力,二次样条就是 s’=0,三次样条就是 s''=0
固定边界 Clamped
首尾两端点的微分值被指定
非节点边界 Not-A-Knot
把端点当做中间点处理,三次函数不做假设,即
三次样条的phthon实现
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)
三次样条函数运行结果如下:
---------------------
作者:宁悦
来源:CSDN
原文:https://blog.csdn.net/deramer1/article/details/79034201
版权声明:本文为作者原创文章,转载请附上博文链接!
内容解析By:CSDN,CNBLOG博客文章一键转载插件