插值與拟合

實際問題中，僅能觀測到函數f(x)在某區間[a,b]上一系列點的值 y i = f ( x i )   , i = 0 , 1 , ⋯   , n y_i=f(x_i)\ ,i=0,1,\cdots,n yi​=f(xi​) ,i=0,1,⋯,n。想得知觀測點之間的函數值時，可用插值或拟合進行近似。

• 插值：用較簡單的、滿足一定條件的函數 φ ( x ) \varphi(x) φ(x)代替 f ( x ) f(x) f(x)，且插值函數滿足 φ ( x i ) = y i   , i = 0 , 1 , ⋯   , n \varphi(x_i)=y_i\ ,i=0,1,\cdots,n φ(xi​)=yi​ ,i=0,1,⋯,n
• 拟合：拟合的近似函數不要求過已知資料點，隻要求在某種意義下在這些點上的總偏差最小

1. 插值方法

1.1 分段線性插值

将每兩個相鄰的節點用直線連起來，如此形成的一條折線就是分段線性插值函數，記為 I n ( x ) I_n(x) In​(x)，滿足 I n ( x i ) = y i I_n(x_i)=y_i In​(xi​)=yi​，且 I n ( x ) I_n(x) In​(x)在每個小區間 [ x i , x i + 1 ] [x_i,x_{i+1}] [xi​,xi+1​]上是線性函數。 I n ( x ) I_n(x) In​(x)可以表示為 I n ( x ) = ∑ i = 0 n y i l i ( x ) I_n(x)=\sum_{i=0}^ny_il_i(x) In​(x)=∑i=0n​yi​li​(x)，其中 l i ( x ) = { x − x i − 1 x i − x i − 1   , x ∈ [ x i − 1 , x i ]   , i ≠ 0 x − x i + 1 x i − x i + 1   , x ∈ [ x i , x i + 1 ]   , i ≠ n 0   , 其 他 l_i(x)=\begin{cases}\cfrac{x-x_{i-1}}{x_i-x_{i-1}}\ ,x\in [x_{i-1},x_i]\ ,i\neq 0 \\ \cfrac{x-x_{i+1}}{x_i-x_{i+1}}\ ,x\in [x_{i},x_{i+1}]\ ,i\neq n \\ 0\ ,其他\end{cases} li​(x)=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​xi​−xi−1​x−xi−1​​ ,x∈[xi−1​,xi​] ,i​=0xi​−xi+1​x−xi+1​​ ,x∈[xi​,xi+1​] ,i​=n0 ,其他​用 I n ( x ) I_n(x) In​(x)計算x點的插值時，隻用到x左右的兩個節點，計算量與節點個數n無關。但n越大，分段越多，插值誤差越小。

1.2 拉格朗日插值多項式

拉格朗日插值的基函數為 l i ( x ) = ∏ j = 0 j ≠ i n x − x j x i − x j   , i = 0 , 1 , ⋯   , n l_i(x)=\prod_{j=0\atop j\neq i}^n \cfrac{x-x_j}{x_i-x_j}\ ,i=0,1,\cdots,n li​(x)=j​=ij=0​∏n​xi​−xj​x−xj​​ ,i=0,1,⋯,n其中 l i ( x ) l_i(x) li​(x)是n次多項式，滿足 l i ( x j ) = { 0   , j ≠ i 1   , j = i l_i(x_j)=\begin{cases}0\ ,j\neq i \\ 1\ ,j=i\end{cases} li​(xj​)={0 ,j​=i1 ,j=i​拉格朗日插值函數為 L n ( x ) = ∑ i = 0 n y i l i ( x ) = ∑ i = 0 n y i ( ∏ j = 0 j ≠ i n x − x j x i − x j ) L_n(x)=\sum_{i=0}^ny_il_i(x)=\sum_{i=0}^ny_i(\prod_{j=0\atop j\neq i}^n \cfrac{x-x_j}{x_i-x_j}) Ln​(x)=i=0∑n​yi​li​(x)=i=0∑n​yi​(j​=ij=0​∏n​xi​−xj​x−xj​​)

1.3 樣條插值

許多工程問題對插值函數的光滑性有較高的要求，要求曲線不僅要連續，而且要有連續的曲率，這就導緻了樣條插值的産生。

1.3.1 樣條函數的概念

數學上将具有一定光滑性的分段多項式稱為樣條函數。具體地說，給定區間[a,b]的一個劃分 Δ : a = x 0 < x 1 < ⋯ < x n − 1 < x n = b \Delta:a=x_0<x_1<\cdots<x_{n-1}<x_n=b Δ:a=x0​<x1​<⋯<xn−1​<xn​=b如果函數S(x)滿足：

• 在每個小區間 [ x i , x i + 1 ] [x_i,x_{i+1}] [xi​,xi+1​]上S(x)是m次多項式
• S(x)在[a,b]上具有m-1階連續導數

則稱S(x)為關于劃分 Δ \Delta Δ的m次樣條函數，其圖形為m次樣條曲線。折線是一次樣條曲線。

1.3.2 三次樣條插值

取插值函數為樣條函數，即為樣條插值。記函數S(x)為f(x)的三次樣條插值函數，且有 S ( x ) = { S i ( x )   , x ∈ [ x i , x i + 1 ]   , i = 0 , 1 , ⋯   , n − 1 } S i ( x ) = a i x 3 + b i x 2 + c i x + d i S(x)=\left\{S_i(x)\ ,x\in [x_i,x_{i+1}]\ ,i=0,1,\cdots,n-1\right\}\\ S_i(x)=a_ix^3+b_ix^2+c_ix+d_i S(x)={Si​(x) ,x∈[xi​,xi+1​] ,i=0,1,⋯,n−1}Si​(x)=ai​x3+bi​x2+ci​x+di​式中 a i , b i , c i , d i a_i,b_i,c_i,d_i ai​,bi​,ci​,di​為待定系數，共4n個。由二階連續可微的性質有 { S i ( x i + 1 ) = S i + 1 ( x i + 1 ) S ’ i ( x i + 1 ) = S ’ i + 1 ( x i + 1 ) S ’ ’ i ( x i + 1 ) = S ’ ’ i + 1 ( x i + 1 ) , i = 0 , 1 , ⋯   , n − 2 \begin{cases}S_i(x_{i+1})=S_{i+1}(x_{i+1})\\ S’_i(x_{i+1})=S’_{i+1}(x_{i+1})\\S’’_i(x_{i+1})=S’’_{i+1}(x_{i+1})\end{cases},i=0,1,\cdots,n-2 ⎩⎪⎨⎪⎧​Si​(xi+1​)=Si+1​(xi+1​)S’i​(xi+1​)=S’i+1​(xi+1​)S’’i​(xi+1​)=S’’i+1​(xi+1​)​,i=0,1,⋯,n−2聯立 S ( x i ) = y i ( i = 0 , 1 , ⋯   , n ) S(x_i)=y_i(i=0,1,\cdots,n) S(xi​)=yi​(i=0,1,⋯,n)共有4n-2個方程，為确定S(x)的4n個待定參數，需再給出兩個邊界條件。

常用的三次樣條函數的邊界條件由3種類型：

1. S ’ ( a ) = y ’ 0 , S ’ ( b ) = y ’ n S’(a)=y’_0,S’(b)=y’_n S’(a)=y’0​,S’(b)=y’n​。由這種邊界條件建立的樣條插值函數稱為f(x)的完備三次樣條插值函數。

   如果f’(x)不知道，可要求S’(x)與f’(x)在端點處近似相等。此時以 x 0 , x 1 , x 2 , x 3 x_0,x_1,x_2,x_3 x0​,x1​,x2​,x3​為節點作一個三次Newton插值多項式 N a ( x ) N_a(x) Na​(x)，以 x n , x n − 1 , x n − 2 , x n − 3 x_n,x_{n-1},x_{n-2},x_{n-3} xn​,xn−1​,xn−2​,xn−3​為節點作一個三次Newton插值多項式 N b ( x ) N_b(x) Nb​(x)，要求 S ’ ( a ) = N a ′ ( a ) , S ’ ( b ) = N ’ b ( b ) S’(a)=N'_a(a),S’(b)=N’_b(b) S’(a)=Na′​(a),S’(b)=N’b​(b)。由這種邊界條件建立的三次樣條稱為f(x)的Lagrange三次樣條插值函數。

2. S ’ ’ ( a ) = y ’ ’ 0 , S ’ ’ ( b ) = y ’ ’ n S’’(a)=y’’_0,S’’(b)=y’’_n S’’(a)=y’’0​,S’’(b)=y’’n​。特别地， y ’ ’ 0 = y ’ ’ n = 0 y’’_0=y’’_n=0 y’’0​=y’’n​=0時，稱為自然邊界條件。
3. S ’ ( a + 0 ) = S ’ ( b − 0 ) , S ’ ’ ( a + 0 ) = S ’ ’ ( b − 0 ) S’(a+0)=S’(b-0),S’’(a+0)=S’’(b-0) S’(a+0)=S’(b−0),S’’(a+0)=S’’(b−0)，此條件稱為周期條件。

1.4 Matlab插值工具箱

1.4.1 一維插值函數

y=interp1(x0,y0,x,’method’)		%其中x0,y0為已知資料點,x為插值點，y為插值點的函數值，method指定插值的方法，預設為線性插值，其值可以為：
	‘nearest’最近項插值	‘linear’線性插值	‘spline’立方樣條插值	‘cubic’立方插值	所有的插值方法要求x0是單調的
           

1.4.2 三次樣條插值

如果三次樣條插值沒有邊界條件，最常用的方法就是采用非扭結（not-a-knot）條件。該條件強迫第一個和第二個三次多項式的三階導數相等，最後一個和倒數第二個三次多項式的三階導數相等。

Matlab中三次樣條插值由如下函數：

y=interp1(x0,y0,x,’spline’)
y=spline(x0,y0,x)
pp=csape(x0,y0)	%提倡使用該函數，csape的傳回值是pp形式，要得到插值點的函數值，需調用函數fnval
	y=fnval(pp,x)，預設邊界條件為Lagrange邊界條件
pp=csape(x0,y0,conds,valconds)	%conds指定插值的邊界條件，當邊界設為二階導數時，valconds給定二階導數的值
           

1.4.3 二維插值

若節點是二維的，插值函數就是二進制函數，即曲面。

• 插值節點為網格節點

z=interp2(x0,y0,z0,x,y,’method’)	%x0,y0分别為m維和n維向量，表示節點；z0為n*m矩陣，表示節點值
%x,y為一維數組，表示插值點；z為矩陣，其行數為y的維數，列數為x的維數，表示得到的插值
pp=csape({x0,y0},z0,conds,valconds),z=fnval(pp,{x,y})	%三次樣條插值
           

• 插值節點為散亂節點

ZI=griddata(x,y,z,XI,YI)	%x,y,z均為n維向量，指明所給資料點的橫、縱、豎坐标
%向量XI,YI為給定網格點的橫坐标和縱坐标，ZI為網格(XI,YI)處的函數值
           

2. 曲線拟合的線性最小二乘法

2.1 線性最小二乘法

已知一組資料，即平面上n個點 ( x i , y i ) (x_i,y_i) (xi​,yi​)，要尋找一個函數f(x)使曲線拟合最好。線性最小二乘法是解決曲線拟合最常用的方法，基本思路是，令 f ( x ) = a 1 r 1 ( x ) + a 2 r 2 ( x ) + ⋯ + a m r m ( x ) f(x)=a_1r_1(x)+a_2r_2(x)+\cdots+a_mr_m(x) f(x)=a1​r1​(x)+a2​r2​(x)+⋯+am​rm​(x)其中 r k ( x ) r_k(x) rk​(x)為事先標明的一組線性無關的函數， a k a_k ak​為待定系數，m<n。拟合準則是使 y i y_i yi​與 f ( x i ) f(x_i) f(xi​)的距離 δ i \delta_i δi​的平方和最小，稱為最小二乘準則。

2.1.1 系數 a k a_k ak​的确定

記 J ( a 1 , ⋯   , a m ) = ∑ i = 1 n δ i 2 = ∑ i = 1 n [ f ( x i ) − y i ] 2 J(a_1,\cdots,a_m)=\sum_{i=1}^n\delta_i^2=\sum_{i=1}^n[f(x_i)-y_i]^2 J(a1​,⋯,am​)=i=1∑n​δi2​=i=1∑n​[f(xi​)−yi​]2要使J達到最小，即要 ∂ J ∂ a j = 0 \cfrac{\partial J}{\partial a_j}=0 ∂aj​∂J​=0，可得到關于 a 1 , ⋯   . a m a_1,\cdots.a_m a1​,⋯.am​的線性方程組，矩陣化後為 R T R A = R T Y \bm{R^TRA=R^TY} RTRA=RTY式中 R = [ r 1 ( x 1 ) ⋯ r m ( x 1 ) ⋮ ⋮ ⋮ r 1 ( x n ) ⋯ r m ( x n ) ] n × m , A = [ a 1 , ⋯   , a m ] T   , Y = [ y 1 , ⋯   , y n ] T \bm{R}=\begin{bmatrix}r_1(x_1) & \cdots & r_m(x_1) \\ \vdots & \vdots & \vdots \\ r_1(x_n) & \cdots & r_m(x_n)\end{bmatrix}_{n\times m},\bm{A}=[a_1,\cdots,a_m]^T\ ,\bm{Y}=[y_1,\cdots,y_n]^T R=⎣⎢⎡​r1​(x1​)⋮r1​(xn​)​⋯⋮⋯​rm​(x1​)⋮rm​(xn​)​⎦⎥⎤​n×m​,A=[a1​,⋯,am​]T ,Y=[y1​,⋯,yn​]T當 { r 1 ( x ) , ⋯   , r m ( x ) } \left\{r_1(x),\cdots,r_m(x)\right\} {r1​(x),⋯,rm​(x)}線性無關時， R \bm{R} R列滿秩， R T R \bm{R^TR} RTR可逆，于是方程組式有唯一解 A = ( R T R ) − 1 R T Y \bm{A=(R^TR)^{-1}R^TY} A=(RTR)−1RTY

2.1.2 函數 r k ( x ) r_k(x) rk​(x)的選取

若通過機理分析，能夠知道y和x之間的函數關系，則 r 1 ( x ) , ⋯   , r m ( x ) r_1(x),\cdots,r_m(x) r1​(x),⋯,rm​(x)容易确定。若無法知道y和x之間的關系，通常可以将資料 ( x i , y i )   , i = 1 , 2 , ⋯   , n (x_i,y_i)\ ,i=1,2,\cdots,n (xi​,yi​) ,i=1,2,⋯,n作圖，直覺地判斷應該用什麼樣的曲線取作拟合。常用的曲線有：

• 直線： y = a 1 x + a 2 y=a_1x+a_2 y=a1​x+a2​
• 多項式： y = a 1 x m + ⋯ + a m x + a m + 1 y=a_1x^m+\cdots+a_mx+a_{m+1} y=a1​xm+⋯+am​x+am+1​
• 雙曲線： y = a 1 x + a 2 y=\cfrac{a_1}{x}+a_2 y=xa1​​+a2​
• 指數曲線： y = a 1 e a 2 x y=a_1e^{a_2x} y=a1​ea2​x，用指數曲線拟合前需作變量代換，化為對 a 1 , a 2 a_1,a_2 a1​,a2​的線性函數

2.2 最小二乘法的Matlab實作

2.2.1 解方程組方法

Matlab中的線性最小二乘标準型為 min ⁡ A ∣ ∣ R A − Y ∣ ∣ 2 2 \mathop{\min}\limits_{A}||\bm{RA-Y}||^2_2 Amin​∣∣RA−Y∣∣22​，求解指令為 A = R ∖ Y \bm{A=R\setminus Y} A=R∖Y。

2.2.2 多項式拟合方法

若取 { r 1 ( x ) , ⋯   , r m + 1 ( x ) } = { 1 , x , ⋯   , x m } \left\{r_1(x),\cdots,r_{m+1}(x)\right\}=\left\{1,x,\cdots,x^m\right\} {r1​(x),⋯,rm+1​(x)}={1,x,⋯,xm}，即用m次多項式拟合給定資料，Matlab中有現成的函數

a=polyfit(x0,y0,m)		%x0,y0為要拟合的資料，m為拟合多項式的次數，a為拟合多項式的系數向量
y=polyval(a,x)		%計算多項式在x出的值
           

3. 最小二乘優化

在非線性規劃問題中，若目标函數由若幹個函數的平方和構成，把極小化這類函數的問題稱為最小二乘優化問題 min ⁡ F ( x ) = ∑ i = 1 m f i 2 ( x ) \min F(x)=\sum_{i=1}^mf_i^2(x) minF(x)=i=1∑m​fi2​(x)在Matlab優化工具箱中，用于求解最小二乘優化問題的函數由lsqlin、lsqcurvefit、lsqnonlin、lsqnonneg。

1. lsqlin函數：求解 min ⁡ x 1 2 ∣ ∣ C ⋅ x − d ∣ ∣ 2 2 s.t. { A ⋅ x ⩽ b A e q ⋅ x = b e q l b ⩽ x ⩽ u b \min_x \cfrac12||\bm{C\cdot x-d}||_2^2 \\ \text{s.t.}\begin{cases} \bm{A\cdot x\leqslant b} \\ Aeq\cdot \bm{x}=beq \\ lb \leqslant\bm{x}\leqslant ub \end{cases} xmin​21​∣∣C⋅x−d∣∣22​s.t.⎩⎪⎨⎪⎧​A⋅x⩽bAeq⋅x=beqlb⩽x⩽ub​式中 C , A , A e q \bm{C,A},Aeq C,A,Aeq為矩陣， d , x , b , b e q , l b , u b \bm{d,x,b},beq,lb,ub d,x,b,beq,lb,ub為向量。

   MATLAB中的函數為

1. lsqcurvefit函數：給定輸入輸出數列xdata，ydata，求參量x，使得 min ⁡ x ∣ ∣ F ( x , x d a t a ) − y d a t a ∣ ∣ 2 2 = ∑ i ( F ( x , x d a t a i ) − y d a t a i ) 2 \min_x||F(x,xdata)-ydata||^2_2=\sum_i(F(x,xdata_i)-ydata_i)^2 xmin​∣∣F(x,xdata)−ydata∣∣22​=i∑​(F(x,xdatai​)−ydatai​)2MATLAB中的函數為

x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)	%fun為定義函數F(x,xdata)的M檔案
           

1. lsqnonlin函數：已知函數向量 F ( x ) = [ f 1 ( x ) , ⋯   , f k ( x ) ] T F(\bm{x})=[f_1(\bm{x}),\cdots,f_k(\bm{x})]^T F(x)=[f1​(x),⋯,fk​(x)]T，求x使得 min ⁡ x ∣ ∣ F ( x ) ∣ ∣ 2 2 \min_x||F(\bm{x})||_2^2 xmin​∣∣F(x)∣∣22​MATLAB中的函數為

x = lsqnonlin(fun,x0,lb,ub,options)		%fun為定義向量函數F(x)的M檔案
           

1. lsqnonneg函數：求解非負的x，使得 min ⁡ x ∣ ∣ C ⋅ x − d ∣ ∣ 2 2 \min_x ||\bm{C\cdot x-d}||_2^2 xmin​∣∣C⋅x−d∣∣22​MATLAB中的函數為

4. 曲線拟合與函數逼近

曲線拟合是選擇一個較簡單函數f(x)在一定準則下最接近一組離散資料。函數逼近是選擇一個較簡單函數f(x)在一定準則下最接近一個較複雜的連續函數y(x),[a,b]。與曲線拟合的最小二乘準則對應，函數逼近常用最小平方逼近，即 J = ∫ a b [ f ( x ) − y ( x ) ] 2 d x J=\int_a^b[f(x)-y(x)]^2dx J=∫ab​[f(x)−y(x)]2dx達到最小。與曲線拟合一樣，選一組函數 { r k ( x ) , k = 1 , ⋯   , m } \left\{r_k(x),k=1,\cdots,m\right\} {rk​(x),k=1,⋯,m}構造f(x)，即令 f ( x ) = a 1 r 1 ( x ) + ⋯ + a m r m ( x ) f(x)=a_1r_1(x)+\cdots+a_mr_m(x) f(x)=a1​r1​(x)+⋯+am​rm​(x)利用極值必要條件可得 [ ( r 1 , r 1 ) ⋯ ( r 1 , r m ) ⋮ ⋮ ⋮ ( r m , r 1 ) ⋯ ( r m , r m ) ] [ a 1 ⋮ a m ] = [ ( y , r 1 ) ⋮ ( y , r m ) ] \left[\begin{array}{ccc} \left(r_{1}, r_{1}\right) & \cdots & \left(r_{1}, r_{m}\right) \\ \vdots & \vdots & \vdots \\ \left(r_{m}, r_{1}\right) & \cdots & \left(r_{m}, r_{m}\right) \end{array}\right]\left[\begin{array}{c} a_{1} \\ \vdots \\ a_{m} \end{array}\right]=\left[\begin{array}{c} \left(y, r_{1}\right) \\ \vdots \\ \left(y, r_{m}\right) \end{array}\right] ⎣⎢⎡​(r1​,r1​)⋮(rm​,r1​)​⋯⋮⋯​(r1​,rm​)⋮(rm​,rm​)​⎦⎥⎤​⎣⎢⎡​a1​⋮am​​⎦⎥⎤​=⎣⎢⎡​(y,r1​)⋮(y,rm​)​⎦⎥⎤​其中 ( g , h ) = ∫ a b g ( x ) h ( x ) d x (g,h)=\int_a^bg(x)h(x)dx (g,h)=∫ab​g(x)h(x)dx，當方程組的系數矩陣非奇異時，有唯一解。