python_matlab_样条插值求未知曲线的函数解析式

Posted 2021-02-14 兵临城下的匹夫

一、问题引入

    对于给出如下的离散的数据点，现在想根据如下的数据点来推测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实现，最后将结果绘制在图上。

  1 import numpy as np
  2 import matplotlib.pyplot as plt
  3 from pylab import mpl
  4 """
  5 三次样条实现：
  6 函数x的自变量为:3, 4.5, 7, 9
  7 因变量为：2.5, 1 2.5, 0.5
  8 """
  9 x = [3, 4.5, 7, 9]
 10 y = [2.5, 1, 2.5, 0.5]
 11 
 12 
 13 """
 14 功能：完后对三次样条函数求解方程参数的输入
 15 参数：要进行三次样条曲线计算的自变量
 16 返回值：方程的参数
 17 """
 18 def calculateEquationParameters(x):
 19 #parameter为二维数组，用来存放参数，sizeOfInterval是用来存放区间的个数
 20 parameter = []
 21 sizeOfInterval=len(x)-1;
 22 i = 1
 23 #首先输入方程两边相邻节点处函数值相等的方程为2n-2个方程
 24 while i < len(x)-1:
 25 data = init(sizeOfInterval*4)
 26 data[(i-1)*4] = x[i]*x[i]*x[i]
 27 data[(i-1)*4+1] = x[i]*x[i]
 28 data[(i-1)*4+2] = x[i]
 29 data[(i-1)*4+3] = 1
 30 data1 =init(sizeOfInterval*4)
 31 data1[i*4] =x[i]*x[i]*x[i]
 32 data1[i*4+1] =x[i]*x[i]
 33 data1[i*4+2] =x[i]
 34 data1[i*4+3] = 1
 35 temp = data[2:]
 36 parameter.append(temp)
 37 temp = data1[2:]
 38 parameter.append(temp)
 39 i += 1
 40 # 输入端点处的函数值。为两个方程, 加上前面的2n - 2个方程，一共2n个方程
 41 data = init(sizeOfInterval * 4 - 2)
 42 data[0] = x[0]
 43 data[1] = 1
 44 parameter.append(data)
 45 data = init(sizeOfInterval * 4)
 46 data[(sizeOfInterval - 1) * 4 ] = x[-1] * x[-1] * x[-1]
 47 data[(sizeOfInterval - 1) * 4 + 1] = x[-1] * x[-1]
 48 data[(sizeOfInterval - 1) * 4 + 2] = x[-1]
 49 data[(sizeOfInterval - 1) * 4 + 3] = 1
 50 temp = data[2:]
 51 parameter.append(temp)
 52 # 端点函数一阶导数值相等为n-1个方程。加上前面的方程为3n-1个方程。
 53 i=1
 54 while i < sizeOfInterval:
 55 data = init(sizeOfInterval * 4)
 56 data[(i - 1) * 4] = 3 * x[i] * x[i]
 57 data[(i - 1) * 4 + 1] = 2 * x[i]
 58 data[(i - 1) * 4 + 2] = 1
 59 data[i * 4] = -3 * x[i] * x[i]
 60 data[i * 4 + 1] = -2 * x[i]
 61 data[i * 4 + 2] = -1
 62 temp = data[2:]
 63 parameter.append(temp)
 64 i += 1
 65 # 端点函数二阶导数值相等为n-1个方程。加上前面的方程为4n-2个方程。且端点处的函数值的二阶导数为零，为两个方程。总共为4n个方程。
 66 i = 1
 67 while i < len(x) - 1:
 68 data = init(sizeOfInterval * 4)
 69 data[(i - 1) * 4] = 6 * x[i]
 70 data[(i - 1) * 4 + 1] = 2
 71 data[i * 4] = -6 * x[i]
 72 data[i * 4 + 1] = -2
 73 temp = data[2:]
 74 parameter.append(temp)
 75 i += 1
 76 return parameter
 77 
 78 
 79 
 80 """
 81 对一个size大小的元组初始化为0
 82 """
 83 def init(size):
 84 j = 0;
 85 data = []
 86 while j < size:
 87 data.append(0)
 88 j += 1
 89 return data
 90 
 91 """
 92 功能：计算样条函数的系数。
 93 参数：parametes为方程的系数，y为要插值函数的因变量。
 94 返回值：三次插值函数的系数。
 95 """
 96 
 97 def solutionOfEquation(parametes,y):
 98 sizeOfInterval = len(x) - 1;
 99 result = init(sizeOfInterval*4-2)
100 i=1
101 while i<sizeOfInterval:
102 result[(i-1)*2]=y[i]
103 result[(i-1)*2+1]=y[i]
104 i+=1
105 result[(sizeOfInterval-1)*2]=y[0]
106 result[(sizeOfInterval-1)*2+1]=y[-1]
107 a = np.array(calculateEquationParameters(x))
108 b = np.array(result)
109 for data_x in b:
110 print(data_x)
111 return np.linalg.solve(a,b)
112 
113 """
114 功能：根据所给参数，计算三次函数的函数值：
115 参数:parameters为二次函数的系数，x为自变量
116 返回值：为函数的因变量
117 """
118 def calculate(paremeters,x):
119 result=[]
120 for data_x in x:
121 result.append(paremeters[0]*data_x*data_x*data_x+paremeters[1]*data_x*data_x+paremeters[2]*data_x+paremeters[3])
122 return result
123 
124 
125 """
126 功能：将函数绘制成图像
127 参数：data_x,data_y为离散的点.new_data_x,new_data_y为由拉格朗日插值函数计算的值。x为函数的预测值。
128 返回值：空
129 """
130 def Draw(data_x,data_y,new_data_x,new_data_y):
131 plt.plot(new_data_x, new_data_y, label="拟合曲线", color="black")
132 plt.scatter(data_x,data_y, label="离散数据",color="red")
133 mpl.rcParams[\'font.sans-serif\'] = [\'SimHei\']
134 mpl.rcParams[\'axes.unicode_minus\'] = False
135 plt.title("三次样条函数")
136 plt.legend(loc="upper left")
137 plt.show()
138 
139 
140 result=solutionOfEquation(calculateEquationParameters(x),y)
141 new_data_x1=np.arange(3, 4.5, 0.1)
142 new_data_y1=calculate([0,0,result[0],result[1]],new_data_x1)
143 new_data_x2=np.arange(4.5, 7, 0.1)
144 new_data_y2=calculate([result[2],result[3],result[4],result[5]],new_data_x2)
145 new_data_x3=np.arange(7, 9.5, 0.1)
146 new_data_y3=calculate([result[6],result[7],result[8],result[9]],new_data_x3)
147 new_data_x=[]
148 new_data_y=[]
149 new_data_x.extend(new_data_x1)
150 new_data_x.extend(new_data_x2)
151 new_data_x.extend(new_data_x3)
152 new_data_y.extend(new_data_y1)
153 new_data_y.extend(new_data_y2)
154 new_data_y.extend(new_data_y3)
155 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实现

  1 import numpy as np
  2 import matplotlib.pyplot as plt
  3 from pylab import mpl
  4 """
  5 三次样条实现：
  6 函数x的自变量为:3,   4.5, 7,    9
  7       因变量为：2.5, 1   2.5,  0.5
  8 """
  9 x = [3, 4.5, 7, 9]
 10 y = [2.5, 1, 2.5, 0.5]
 11  
 12  
 13 """
 14 功能：完后对三次样条函数求解方程参数的输入
 15 参数：要进行三次样条曲线计算的自变量
 16 返回值：方程的参数
 17 """
 18 def calculateEquationParameters(x):
 19     #parameter为二维数组，用来存放参数，sizeOfInterval是用来存放区间的个数
 20     parameter = []
 21     sizeOfInterval=len(x)-1;
 22     i = 1
 23     #首先输入方程两边相邻节点处函数值相等的方程为2n-2个方程
 24     while i < len(x)-1:
 25         data = init(sizeOfInterval*4)
 26         data[(i-1)*4] = x[i]*x[i]*x[i]
 27         data[(i-1)*4+1] = x[i]*x[i]
 28         data[(i-1)*4+2] = x[i]
 29         data[(i-1)*4+3] = 1
 30         data1 =init(sizeOfInterval*4)
 31         data1[i*4] =x[i]*x[i]*x[i]
 32         data1[i*4+1] =x[i]*x[i]
 33         data1[i*4+2] =x[i]
 34         data1[i*4+3] = 1
 35         temp = data[2:]
 36         parameter.append(temp)
 37         temp = data1[2:]
 38         parameter.append(temp)
 39         i += 1
 40    # 输入端点处的函数值。为两个方程, 加上前面的2n - 2个方程，一共2n个方程
 41     data = init(sizeOfInterval * 4 - 2)
 42     data[0] = x[0]
 43     data[1] = 1
 44     parameter.append(data)
 45     data = init(sizeOfInterval * 4)
 46     data[(sizeOfInterval - 1) * 4 ] = x[-1] * x[-1] * x[-1]
 47     data[(sizeOfInterval - 1) * 4 + 1] = x[-1] * x[-1]
 48     data[(sizeOfInterval - 1) * 4 + 2] = x[-1]
 49     data[(sizeOfInterval - 1) * 4 + 3] = 1
 50     temp = data[2:]
 51     parameter.append(temp)
 52     # 端点函数一阶导数值相等为n-1个方程。加上前面的方程为3n-1个方程。
 53     i=1
 54     while i < sizeOfInterval:
 55         data = init(sizeOfInterval * 4)
 56         data[(i - 1) * 4] = 3 * x[i] * x[i]
 57         data[(i - 1) * 4 + 1] = 2 * x[i]
 58         data[(i - 1) * 4 + 2] = 1
 59         data[i * 4] = -3 * x[i] * x[i]
 60         data[i * 4 + 1] = -2 * x[i]
 61         data[i * 4 + 2] = -1
 62         temp = data[2:]
 63         parameter.append(temp)
 64         i += 1
 65     # 端点函数二阶导数值相等为n-1个方程。加上前面的方程为4n-2个方程。且端点处的函数值的二阶导数为零，为两个方程。总共为4n个方程。
 66     i = 1
 67     while i < len(x) - 1:
 68         data = init(sizeOfInterval * 4)
 69         data[(i - 1) * 4] = 6 * x[i]
 70         data[(i - 1) * 4 + 1] = 2
 71         data[i * 4] = -6 * x[i]
 72         data[i * 4 + 1] = -2
 73         temp = data[2:]
 74         parameter.append(temp)
 75         i += 1
 76     return parameter
 77  
 78  
 79  
 80 """
 81 对一个size大小的元组初始化为0
 82 """
 83 def init(size):
 84     j = 0;
 85     data = []
 86     while j < size:
 87         data.append(0)
 88         j += 1
 89     return data
 90  
 91 """
 92 功能：计算样条函数的系数。
 93 参数：parametes为方程的系数，y为要插值函数的因变量。
 94 返回值：三次插值函数的系数。
 95 """
 96  
 97 def solutionOfEquation(parametes,y):
 98     sizeOfInterval = len(x) - 1;
 99     result = init(sizeOfInterval*4-2)
100     i=1
101     while i<sizeOfInterval:
102         result[(i-1)*2]=y[i]
103         result[(i-1)*2+1]=y[i]
104         i+=1
105     result[(sizeOfInterval-1)*2]=y[0]
106     result[(sizeOfInterval-1)*2+1]=y[-1]
107     a = np.array(calculateEquationParameters(x))
108     b = np.array(result)
109     for data_x in b:
110         print(data_x)
111     return np.linalg.solve(a,b)
112  
113 """
114 功能：根据所给参数，计算三次函数的函数值：
115 参数:parameters为二次函数的系数，x为自变量
116 返回值：为函数的因变量
117 """
118 def calculate(paremeters,x):
119     result=[]
120     for data_x in x:
121         result.append(paremeters[0]*data_x*data_x*data_x+paremeters[1]*data_x*data_x+paremeters[2]*data_x+paremeters[3])
122     return  result
123  
124  
125 """
126 功能：将函数绘制成图像
127 参数：data_x,data_y为离散的点.new_data_x,new_data_y为由拉格朗日插值函数计算的值。x为函数的预测值。
128 返回值：空
129 """
130 def  Draw(data_x,data_y,new_data_x,new_data_y):
131         plt.plot(new_data_x, new_data_y, label="拟合曲线", color="black")
132         plt.scatter(data_x,data_y, label="离散数据",color="red")
133         mpl.rcParams[\'font.sans-serif\'] = [\'SimHei\']
134         mpl.rcParams[\'axes.unicode_minus\'] = False
135         plt.title("三次样条函数")
136         plt.legend(loc="upper left")
137         plt.show()
138  
139  
140 result=solutionOfEquation(calculateEquationParameters(x),y)
141 new_data_x1=np.arange(3, 4.5, 0.1)
142 new_data_y1=calculate([0,0,result[0],result[1]],new_data_x1)
143 new_data_x2=np.arange(4.5, 7, 0.1)
144 new_data_y2=calculate([result[2],result[3],result[4],result[5]],new_data_x2)
145 new_data_x3=np.arange(7, 9.5, 0.1)
146 new_data_y3=calculate([result[6],result[7],result[8],result[9]],new_data_x3)
147 new_data_x=[]
148 new_data_y=[]
149 new_data_x.extend(new_data_x1)
150 new_data_x.extend(new_data_x2)
151 new_data_x.extend(new_data_x3)
152 new_data_y.extend(new_data_y1)
153 new_data_y.extend(new_data_y2)
154 new_data_y.extend(new_data_y3)
155 Draw(x,y,new_data_x,new_data_y)

三次插值

   三次样条函数运行结果如下：
 

 

原文转自：https://blog.csdn.net/deramer1/article/details/79034201

参考博客：https://www.cnblogs.com/ondaytobewhoyouwant/p/8989497.html

                  https://blog.csdn.net/flyingleo1981/article/details/53008931    c语言实现 

                  https://blog.csdn.net/guanmao4322/article/details/85054189    原理及matlab实现

                  https://blog.csdn.net/zb1165048017/article/details/48311603    原理及matlab实现(有手稿哟)

                  https://blog.csdn.net/u012856866/article/details/23952819       c++实现