一、問題引入
對於給出如下的離散的數據點，現在想根據如下的數據點來推測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++實現