|
一、问题引入 对于给出如下的离散的数据点，现在想根据如下的数据点来推测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 三次样条的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 版权声明：本文为博主原创文章，转载请附上博文链接！ |
|
|
