數(shù)模第09章插值與擬合

1、第九章插值與擬合插值：求過已知有限個數(shù)據(jù)點的近似函數(shù)。擬合：已知有限個數(shù)據(jù)點，求近似函數(shù)，不要求過已知數(shù)據(jù)點，只要求在某種意義下它在這些點上的總偏差最小。插值和擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個函數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。而面對一個實際問題，究竟應(yīng)該用插值還是擬合，有時容易確定，有時則并不明顯。1 插值方法下面介紹幾種基本的、常用的插值：拉格朗日多項式插值、牛頓插值、分段線性插值、Hermite 插值和三次樣條插值。拉格朗日多項式插值插值多項式用多項式作為研究插值的工具，稱為代數(shù)插值。其基本問題是：已知函數(shù) f ( x) 在區(qū)間a,b 上n 1 個不同點至多n

2、次多項式n 處的函數(shù)值 yi f ( xi ) (i 0,1,L, n) ，求一個n ( x) a0 a1 x L an xn（1）使其在給定點處與 f ( x) 同值，即滿足插值條件n ( xi ) f ( x)（2）n ( x) 稱為插值多項式， xi (i 0,1,L, n) 稱為插值節(jié)點，簡稱節(jié)點，a,b 稱為插值區(qū)間。從幾何上看， n 次多項式插值就是過n 1 個點( xi , f ( xi ) (i 0,1,L, n) ，作一條多項式曲線 y n ( x) 近似曲線 y f ( x) 。n 次多項式（1）有n 1 個待定系數(shù)，由插值條件（2）恰好給出n 1 個方程a a x a x

3、 L a y2xn01 02 0n 00a a x a x L a y2xn01 12 1n 11（3）LLLLLLLLLLLLaa x a x L a y2xn 01 n2 nn nn記此方程組的系數(shù)矩陣為 A ，則n0n111det( A) LLLLLLLn n1是范德蒙特(Vandermonde)行列式。當(dāng)n 互不相同時，此行列式值不為零。因此方程組（3）有唯一解。這表明，只要n 1 個節(jié)點互不相同，滿足插值要求（2）的插值多項式（1）是唯一的。插值多項式與函數(shù)之間的差Rn ( x) f ( x) n ( x)-175-稱為截斷誤差，又稱為插值余項。當(dāng) f ( x) 充分光滑時，f (

4、n 1) ( )n1( x), (a, b)Rn ( x) f ( x) Ln ( x) nj ) 。(n 1)!其中n1(j 01.1.2拉格朗日插值多項式實際上比較方便的作法不是解方程（3）求待定系數(shù)，而是先構(gòu)造一組基函數(shù)(li ( x) ()x x jn(i 0,1,L,n),x xj0 jiijli ( x) 是n 次多項式，滿足j ij il ( x ) 0ij1令x xnnnL ( x) y l ( x) yj （4）i ni ix xij i 0i 0j 0j i上式稱為n 次 Lagrange 插值多項式，由方程（3）解的唯一性， n 1 個節(jié)點的n 次Lagrange 插值多

5、項式存在唯一。1.1.3用作 Lagrange 插值中沒有現(xiàn)成的Lagrange 插值函數(shù)，必須編寫一個M 文件實現(xiàn)Lagrange 插值。設(shè)n 個節(jié)點數(shù)據(jù)以數(shù)組 x0, y0 輸入（注意 Mat 的數(shù)組下標(biāo)從 1 開始），m 個插值點以數(shù)組 x 輸入，輸出數(shù)組 y 為m 個插值。編寫一個名為 lagrange.m 的 M 文件：function y=lagrange(x0,y0,x); n=length(x0);m=length(x);for i=1:m z=x(i); s=0.0;for k=1:n p=1.0;for j=1:n if j=kp=p*(z-x0(j)/(x0(k)-x0(

6、j); endend s=p*y0(k)+s;end y(i)=s;end-176-12.牛頓（Newton）插值在導(dǎo)出 Newton 公式前，先介紹公式表示中所需要用到的差商、差分的概念及性質(zhì)。1.2.1 差商定義 設(shè)有函數(shù) f (1, x2 ,L為一系列互不相等的點，稱f ( xi ) f ( x j ) (i j) 為 f ( x) 關(guān)于點 x , x 一階差商（也稱均差）記為 f x , x ，即ijijx xijf ( xi ) f ( x j )f x , x ijx xij稱一階差商的差商f xi , x j f x j , xk xi xkk 的二階差商，記為 f 為 f (

7、x) 關(guān)于點f k 。一般地，稱k 1 f k x0 xk為 f ( x) 關(guān)于點k 的k 階差商，記為k 1 f f k f kx x0k容易證明，差商具有下述性質(zhì)：f xi , x j f x j , xi k f j f f k 1.2.2Newton 插值公式線性插值公式可表成1( x) f (0 ) f x0 , x1 稱為一次 Newton 插值多項式。一般地，由各階差商的定義，依次f ( x) f (f x, x0 f 0 ) f x, x0 x1 ) f 2 ( x x2 ) f 11 f f 1, x2 LLn1 n ( x xn ) f 1),L, (f f n 將以上各式

8、分別乘以1,(n1) ，然后相加并消去兩邊相等的部分，即得f ( x) f ( ( (0 ) f x0 , x1 Ln1) f n ) f n 1,L, xn 記-177-Nn ( x) f ( (Rn (0 ) f x0 , x1 Ln1) f xn ) f 1,L, xn n 1,L, xn 0 )( n1( x) f 顯然， Nn ( x) 是至多n 次的多項式，且滿足插值條件，因而它是 f ( x) 的n 次插值多項式。這種形式的插值多項式稱為 Newton 插值多項式。 Rn ( x) 稱為 Newton 插值余項。Newton 插值的優(yōu)點是：每增加一個節(jié)點，插值多項式只增加一項，即

9、Nn 1 ( x) Nn (0 )L( x xn ) f n 1 因而便于遞推運(yùn)算。而且 Newton 插值的計算量小于 Lagrange 插值。由插值多項式的唯一性可知，Newton 插值余項與 Lagrange 余項也是相等的，即Rn ( x) n1( x) f 1,L, xn f ( n 1) ( ) n1( x) (a, b)(n 1)!由此差商與導(dǎo)數(shù)的關(guān)系f (n) ( )n f n!其中 (, ), mxi 。01.2.3差分當(dāng)節(jié)點等距時，即相鄰兩個節(jié)點之差（稱為步長）為常數(shù)，Newton 插值公式的形式會更簡單。此時關(guān)于節(jié)點間函數(shù)的平均變化率（差商）可用函數(shù)值之差（差分）來表示。

10、定義 設(shè)有等距節(jié)點 xk x0 k稱相鄰兩個節(jié)點 xk , xk 1 處的函數(shù)值的增量) ，步長h 為常數(shù)， fk f ( xk ) 。fk (k 0,1,L, n 1) 為函數(shù) f ( x) 在fk 1 點 xk 處以h 為步長的一階差分，記為fk ，即(k 0,1,L, n)k類似地，定義差分的差分為高階差分。如二階差分為2一般地， m 階差分為m(k 0,1,L, n 2)k(k 2,3,L) ，k上面定義的各階差分又稱為向前差分。常用的差分還有兩種：k 1稱為 f ( x) 在 xk 處以h 為步長的向后差分；f f x h f x h 2 2 kkk稱為 f ( x) 在 xk 處以

11、h 為步長的中心差分。一般地， m 階向后差分與m 階中心差分公式為-178-m mk 1k 12k 2差分具有以下性質(zhì)：（i）各階差分均可表成函數(shù)值的線性組合，例如j m mmfk (1) j fk m jj 0j m mmfk (1) j fk jj 0（ii）各種差分之間可以互化。向后差分與中心差分化成向前差分的公式如下：m f m fkk m m f m fkk m21.2.4等距節(jié)點插值公式如果插值節(jié)點是等距的，則插值公式可用差分表示。設(shè)已知節(jié)點xk x0 kh (k 0,1,2,L, n) ，則有Nn ( x) f ( x0 ) f L f x0 )n (n1 )n1 )fn f

12、f0 0 ( x x0 ) L 0 (n!hnh若令 x x0 th ，則上式又可變形為L t(t 1)L(t n 1) n fN ( x th) f tfn0000n!上式稱為 Newton 向前插值公式。分段線性插值插值多項式的振蕩用 Lagrange 插值多項式 Ln ( x) 近似 f ( x)(a x b) ，雖然隨著節(jié)點個數(shù)的增加， Ln ( x) 的次數(shù)n 變大，多數(shù)情況下誤差| Rn ( x) |會變小。但是n 增大時， Ln ( x) 的光滑性變壞，有時會出現(xiàn)很大的振蕩。理論上，當(dāng)n ，在a,b 內(nèi)并不能保證 Ln ( x) 處處收斂于 f ( x) 。Runge 給出了一個

13、有名的例子：1f ( x) ,x 5,51 x2對于較大的| x | ，隨著n 的增大，Ln ( x) 振蕩越來越大，事實上可以證明，僅當(dāng)| x | 3.63時，才有l(wèi)im Ln ( x) f ( x) ，而在此區(qū)間外， Ln ( x) 是發(fā)散的。n高次插值多項式的這些缺陷，促使人們轉(zhuǎn)而尋求簡單的低次多項式插值。1.3.2分段線性插值簡單地說，將每兩個相鄰的節(jié)點用直線連起來，如此形成的一條折線就是分段線性插值函數(shù)，記作 In ( x) ，它滿足 I，且 In ( x) 在每個小區(qū)間xi , xi 1 上是線性-179-函數(shù)(i 0,1,L, n) 。In ( x) 可以表示為In ( x) y

14、ili ( x)ni 0 x xi1 (i 0時舍去）, xixii1 x xl ( x) i1（i n時舍去）， xii1xii10，其它In ( x) 有良好的收斂性，即對于 x a,b 有，lim In ( x) f ( x) 。n用 In ( x) 計算 x 點的插值時，只用到 x 左右的兩個節(jié)點，計算量與節(jié)點個數(shù)n 無關(guān)。但n 越大，分段越多，插值誤差越小。實際上用函數(shù)表作插值計算時，分段線性插值就足夠了，如數(shù)學(xué)、物理中用的特殊函數(shù)表，數(shù)理統(tǒng)計中用的概率分布表等。1.3.3用用1。y=實現(xiàn)分段線性插值實現(xiàn)分段線性插值不需要編制函數(shù)程序，中有現(xiàn)成的一維插值函數(shù)1(x0,y0,x,met

15、hod)method 指定插值的方法，默認(rèn)為線性插值。其值可為：nearest linear spline cubic最近項插值線性插值逐段 3 次樣條插值保凹凸性 3 次插值。所有的插值方法要求 x0 是單調(diào)的。當(dāng) x0 為等距時可以用快速插值法，使用快速插值法的格式為*nearest、*linear、*spline、*cubic。埃爾米特(Hermite)插值Hermite 插值多項式如果對插值函數(shù)，不僅要求它在節(jié)點處與函數(shù)同值，而且要求它與函數(shù)有相同的一階、二階甚至更高階的導(dǎo)數(shù)值，這就是 Hermite 插值問題。本節(jié)主要函數(shù)與函數(shù)的值及一階導(dǎo)數(shù)值均相等的 Hermite 插值。在節(jié)點處

16、插值設(shè)已知函數(shù) y f ( x) 在n 1 個互異節(jié)點n 上的函數(shù)值 yi f ( xi )(i 0,1,L, n) 和導(dǎo)數(shù)值 yi f ( xi ) (i 0,1,L, n) ，要求一個至多2n 1次的多項式H ( x) ，使得H ( xi ) yiH ( xi ) yi(i 0,1,L, n)滿足上述條件的多項式 H ( x) 稱為 Hermite 插值多項式。Hermite 插值多項式為-180-nH ( x) hi ( xi x)(2ai yi yi ) yi i 02 x x j nn1 ， ai 其中hi。x xj 0 xi x jj 0 j iijj i1.4.2用實現(xiàn) Herm

17、ite 插值中沒有現(xiàn)成的 Hermite 插值函數(shù)，必須編寫一個 M 文件實現(xiàn)插值。設(shè)n 個節(jié)點的數(shù)據(jù)以數(shù)組 x0 （已知點的橫坐標(biāo)）， y0 （函數(shù)值）， y1 （導(dǎo)數(shù)值）輸入（注意 M at 的數(shù)組下標(biāo)從 1 開始），m 個插值點以數(shù)組 x 輸入，輸出數(shù)組 y 為m個插值。編寫一個名為 hermite.m 的 M 文件：function y=hermite(x0,y0,y1,x);n=length(x0);m=length(x); for k=1:myy=0.0;for i=1:n h=1.0;a=0.0;for j=1:n if j=ih=h*(x(k)-x0(j)/(x0(i)-x0(

18、j)2; a=1/(x0(i)-x0(j)+a;end endyy=yy+h*(x0(i)-x(k)*(2*a*y0(i)-y1(i)+y0(i); endy(k)=yy; end1.5樣條插值許多工程技術(shù)中計算問題對插值函數(shù)的光滑性有較高要求，如飛機(jī)的機(jī)翼外形，內(nèi)燃機(jī)的進(jìn)、排氣門的凸輪曲線，都要求曲線具有較高的光滑程度，不僅要連續(xù)，而且要有連續(xù)的曲率，這就導(dǎo)致了樣條插值的產(chǎn)生。1.5.1 樣條函數(shù)的概念所謂樣條（Spline）本來是工程設(shè)計中使用的一種繪圖工具，它是富有彈性的細(xì)木條或細(xì)金屬條。繪圖員利用它把一些已知點連接成一條光滑曲線（稱為樣條曲線），并使連接點處有連續(xù)的曲率。數(shù)學(xué)上將具有一

19、定光滑性的分段多項式稱為樣條函數(shù)。具體地說，給定區(qū)間a,b的一個分劃 :a 如果函數(shù) s(x) 滿足：（i）在每個小區(qū)間xn 1 xn b 1) 上 s(x) 是 k 次多項式；（ii） s(x) 在a,b 上具有 k 1 階連續(xù)導(dǎo)數(shù)。則稱 s(x) 為關(guān)于分劃 的k 次樣條函數(shù)，其圖形稱為k 次樣條曲線。n 稱為n1 稱為內(nèi)節(jié)點， x0 , xn 稱為邊界點，這類樣條函數(shù)的全體記做樣條節(jié)點，-181-SP (, k ) ，稱為 k 次樣條函數(shù)空間。顯然，折線是一次樣條曲線。若 s(x) SP (, k ) ，則 s(x) 是關(guān)于分劃 的 k 次多項式樣條函數(shù)。k 次多項式樣條函數(shù)的一般形式為

20、n1 j xiksk (x) i (x xk) ，j i!k!i0j1其中i (i 0,1,L, k ) 和 j ( j 1,2,L, n 1) 均為任意常數(shù)，而( x j，（ j 1,2,L, n 1 ）(x x)kj x x0,j在實際中最常用的是 k 2 和 3 的情況，即為二次樣條函數(shù)和三次樣條函數(shù)。二次樣條函數(shù)：對于a, b 上的分劃 : a b ，則nx2 jn1s (x) x (x x)2 S (,2) ，2（5）j P2012!2!j1( x j其中(x x) 2，（ j 1,2,L, n 1 ）。j x x0,j三次樣條函數(shù)：對于a, b 上的分劃 : a b ，則n3 jn

21、1s (x) (x x )3 S (,3) ，（6）j P3013!j1( x j其中(x x) 3，（ j 1,2,L, n 1 ）。j x x0,j利用樣條函數(shù)進(jìn)行插值，即取插值函數(shù)為樣條函數(shù)，稱為樣條插值。例如分段線性插值是一次樣條插值。下面介紹二次、三次樣條插值。1.5.2二次樣條函數(shù)插值首先，注意到 s2 (x) SP (,2) 中含有 n 2 個特定常數(shù)，故應(yīng)需要 n 2 個插值條件，因此，二次樣條插值問題可分為兩類：問題（1）：已知插值節(jié)點 xi 和相應(yīng)的函數(shù)值 yi (i 0,1,L, n) 以及端點 x0 （或 xn ）處的導(dǎo)數(shù)值 y0 （或 yn ），求 s2 (x) SP

22、 (,2) 使得s2 (xi ) yi (i 0,1,2,L, n)s (x ) y (或s (x ) y )（7） 200nnn問題（2）：已知插值節(jié)點 xi 和相應(yīng)的導(dǎo)數(shù)值 yi (i 0,1,2,L, n) 以及端點 x0 （或 xn ）處的函數(shù)值 y0 （或 yn ），求 s2 (x) SP (,2) 使得s2 (xi ) yi (i 0,1,2,L, n)s（8）(x ) y (或s (x ) y ) 2002nn-182-事實上，可以證明這兩類插值問題都是唯一可解的。對于問題（1），由條件（7）s(x ) x 1 x2 y2001 02 0021s(x ) x yx2 2101 1

23、2 112j111j is (x ) x x (x x ) y j( j 2,3,L, n)222j01 j2ji22i1s2 (x0 ) 1 2 x0 y0引入記號 X ( , , , ,L, )T 為未知向量，C ( y , y ,L, y , y ) 為已n1012101n0知向量。11 x2)2 x00L LL1002122x1x001A 1 x )20M21MMMM10n11x000L (0 ,1,2 , 1,L, n1)T于是，問題轉(zhuǎn)化為求方程組 AX C 的解 X即到二次樣條函數(shù) s2 (x) 的表達(dá)式。對于問題（2）的情況類似。1.5.3三次樣條函數(shù)插值，由于 s3 (x) S

24、P (,3) 中含有 n 3 個待定系數(shù)，故應(yīng)需要 n 3 個插值條件，已知插值節(jié)點 xi 和相應(yīng)的函數(shù)值 f (xi ) yi (i 0,1,2,L, n) ，這里提供了 n 1個條件，還需要 2 個邊界條件。常用的三次樣條函數(shù)的邊界條件有 3 種類型：（i）s3 (a) y0 , s3 (b) yn 。由這種邊界條件建立的樣條插值函數(shù)稱為 f ( x) 的完備三次樣條插值函數(shù)。特別地， y0 yn 0 時，樣條曲線在端點處呈水平狀態(tài)。如果 f ( x) 不知道，可以要求 s3 (x) 與 f ( x) 在端點處近似相等。這時以2 , x3 為節(jié)點作一個三次 Newton 插值多項式 Na

25、( x) ，以個三次 Newton 插值多項式 Nb ( x) ，要求s(a) N a (a), s(b) N b (b)n 2 , xn 3 作一由這種邊界條件建立的三次樣條稱為 f ( x) 的 Lagrange 三次樣條插值函數(shù)。（ii） s3 (a) y0 , s3 (b) y3 。特別地 yn yn 0 時，稱為自然邊界條件。-183-（iii） s3 (a 0) s3 (b 0), s3 (a 0) s3 (b 0) ，(這里要求 s3 (a 0) s3 (b 0) )此條件稱為周期條件。1.5.4 三次樣條插值在中的實現(xiàn)在中數(shù)據(jù)點稱之為斷點。如果三次樣條插值沒有邊界條件，最常用的

26、方法，就是采用非扭結(jié)（not-a-knot）條件。這個條件強(qiáng)迫第 1 個和第 2 個三次多項式的三階導(dǎo)數(shù)相等。對最后一個和倒數(shù)第 2 個三次多項式也做同樣地處理。中三次樣條插值也有現(xiàn)成的函數(shù)：y=1(x0,y0,x,spline)；y=spline(x0,y0,x)；pp=cs(x0,y0,conds)，y=ppval(pp,x)。其中 x0,y0 是已知數(shù)據(jù)點，x 是插值點，y 是插值點的函數(shù)值。對于三次樣條插值，提倡使用函數(shù) cs值點的函數(shù)值，必須調(diào)用函數(shù) ppval。，cs的返回值是 pp 形式，要求插pp=cs pp=cs(x0,y0)：使用默認(rèn)的邊界條件，即 Lagrange 邊界條

27、件。(x0,y0,conds)中的 conds 指定插值的邊界條件，其值可為：complete not-a-knot periodic second variational邊界為一階導(dǎo)數(shù)，即默認(rèn)的邊界條件非扭結(jié)條件周期條件邊界為二階導(dǎo)數(shù)，二階導(dǎo)數(shù)的值0, 0。設(shè)置邊界的二階導(dǎo)數(shù)值為0,0。對于一些特殊的邊界條件，可以通過 conds 的一個1 2 矩陣來表示，conds 元素的取值為 1，2。此時，使用命令pp=cs(x0,y0_ext,conds)其中 y0_ext=left, y0, right，這里 left 表示左邊界的取值，right 表示右邊界的取值。conds(i)=j 的含義是

28、給定端點i 的 j 階導(dǎo)數(shù)，即 conds 的第一個元素表示左邊界的條件，第二個元素表示右邊界的條件，conds=2,1表示左邊界是二階導(dǎo)數(shù)，右邊界是一階導(dǎo)數(shù)，對應(yīng)的值由 left 和 right 給出。詳細(xì)情況請使用幫助 help cs例 1機(jī)床加工。待加工零件的外形根據(jù)工藝要求由一組數(shù)據(jù)( x, y) 給出（在平面情況下），用程控銑床加工時每一刀只能沿 x 方向和 y 方向走非常小的一步，這就需要從已知數(shù)據(jù)得到加工所要求的步長很小的( x, y) 坐標(biāo)。表 1 中給出的 x, y 數(shù)據(jù)位于機(jī)翼斷面的下輪廓線上，假設(shè)需要得到 x 坐標(biāo)每改變0.1 時的 y 坐標(biāo)。試完成加工所需數(shù)據(jù)，畫出曲線

29、，并求出 x 0 處的曲線斜率和13 x 15 范圍內(nèi) y 的最小值。表 1x0357911121314 15y0 1.2 1.7 2.0 2.1 2.0 1.8 1.21.0 1.6要求用Lagrange、分段線性和三次樣條三種插值方法計算。解 編寫以下程序：clc,clearx0=0 3-184-57911121314 15;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 x=0:0.1:15;1.0 1.6;y1=lagrange(x0,y0,x); %調(diào)用前面編寫的Lagrange插值函數(shù)y2= y3= pp1=cs pp2=cs1(x0,y0,x);1(x0,y0,

30、x,spline);(x0,y0); y4=ppval(pp1,x);(x0,y0,second); y5=ppval(pp2,x);fprf(比較一下不同插值方法和邊界條件的結(jié)果:n)fprf(xy1y2y3y4y5n)xianshi=x,y1,y2,y3,y4,y5;fprf(%ft%ft%ft%ft%ft%fn,xianshi) subplot(2,2,1), plot(x0,y0,+,x,y1), title(Lagrange)subplot(2,2,2), plot(x0,y0,+,x,y2), title(Piecewise linear) subplot(2,2,3), plot

31、(x0,y0,+,x,y3), title(Spline1)subplot(2,2,4), plot(x0,y0,+,x,y4), title(Spline2)p1),x0(1) %求x=0處的導(dǎo)數(shù)dyx0=ppval(fndytemp=y3(131:151);index=find(ytemp=min(ytemp); xymin=x(130+index),ytemp(index)計算結(jié)果略?？梢钥闯?，拉格朗日插值的結(jié)果根本不能應(yīng)用，分段線性插值的光滑性較差（特別是在 x 14 附近彎曲處），建議選用三次樣條插值的結(jié)果。B 樣條函數(shù)插值方法磨光函數(shù)實際中的許多問題，往往是既要求近似函數(shù)（曲線或曲

32、面）有足夠的光滑性，又要求與實際函數(shù)有相同的凹凸性，一般插值函數(shù)和樣條函數(shù)都不具有這種性質(zhì)。如果對于一個特殊函數(shù)進(jìn)行磨光處理生成磨光函數(shù)（多項式），則用磨光函數(shù)構(gòu)造出樣條函數(shù)作為插值函數(shù)，既有足夠的光滑性，而且也具有較好的保凹凸性，因此磨光函數(shù)在一維插值（曲線）和二維插值（曲面）問題中有著廣泛的應(yīng)用。由積分理論可知，對于可積函數(shù)通過積分會提高函數(shù)的光滑度，因此，用積分方法對函數(shù)進(jìn)行磨光處理。定義 若 f (x) 為可積函數(shù)，對于h 0 ，則稱積分可以利x h1h f(x) 2 f (t)dt1,hhx2為 f (x) 的一次磨光函數(shù)， h 稱為磨光寬度。 同樣的，可以定義 f (x) 的k 次

33、磨光函數(shù)為x h1h f(x) (t)dt （ k 1 ）2 fh k 1,h2k ,hx事實上，磨光函數(shù) fk ,h (x) 比 f (x) 的光滑程度要高，且當(dāng)磨光寬度h 很少時 fk ,h (x)很接近于 f (x) 。1.6.2等距B 樣條函數(shù)對于任意的函數(shù) f (x) ，定義其步長為 1 的中心差分算子 如下：-185-f (x) f (x 1 ) f (x 1 )22在此取 f (x) x0 ，則01 02 是一個方波函數(shù)（如圖 1），記 (x) x 。并取 h 1，對 (x) 進(jìn)行一次磨光000得x 1 1 0 1 0 x 1 (x) 2(t)dt t t dt211012 2

34、x2x2 x1xxt dt t 0dt ( 1)0 x1顯然1(x) 是連續(xù)的（如圖 1）。圖 10 (x) 和1(x) 的圖形類似地到 k 次磨光函數(shù)為kk 1C jk 1 (x) (x j k 1k!j 1)k2j0實際上，可以證明： k (x) 是分段k 次多項式，且具有k 1階連續(xù)導(dǎo)數(shù)，其k 階k 1導(dǎo)數(shù)有 k 2 個間斷點，記為 x j j 對應(yīng)于分劃 : （ j 0,1,2,L, k 1 ）。從而可知k (x) 是2k1 的 k 次多項式樣條函數(shù)，稱之為基本樣k 1條函數(shù)，簡稱為 k 次 B 樣條。由于樣條節(jié)點為 x j j 等距的，故k (x) 又稱為 k 次等距 B 樣條函數(shù)。

35、（ j 0,1,2,L, k 1 ）是2對于任意函數(shù) f (x) 的 k 次磨光函數(shù)，由歸納法可以得到： x t f(x) 1hhf (t)dt （ x t x ）k 1 k ,hhh22 x t 特別地，當(dāng) f (x) 1 時，有 1dt 1，從而(x)dx 1 ，且當(dāng)k 1 khhk 1時有遞推關(guān)系-186- (x) 1 x k 1 x 1 k 1 x x 1 k k 1 2 k 1 2 k221.6.3一維等距 B 樣條函數(shù)插值等距 B 樣條函數(shù)與通常的樣條有如下的關(guān)系：定理 設(shè)有區(qū)間a, b 的均勻分劃 : x x jh（ j 0,1,2,L, n ），h b a ，j0n則對任意k

36、次樣條函數(shù) sk (x) SP (, k ) 都可以表示為 B 樣條函數(shù)族k 1 jn1 x xk 0 j h2 j k的線性組合。一組點 (x j , y j ) ，其中 x j x0 jh根據(jù)定理 ，如果已 知曲線上 （ h 0, j 0,1,2,L, n ），則可以構(gòu)造出一條樣條磨光曲線（即為 B 樣條函數(shù)族的線性組合）n1 x x0j js (x) c k khj k其中cj （ j k,k 1,L, n 1 ）為待定常數(shù)。用它來近曲線，既有較好的精度，又有良好的保凸性。實際中，最常用的是k 3 的情況，即一般形式為 x x0n1j js (x) c 33hj1其中 n 3 個待定系數(shù)

37、cj （ j 1,0,L, n 1）可以由插值條件確定。對于插值條件s3 (x j ) y j ( j 0,1,2,Ln)s (x ) y ( j 0, n)3jj有sn(x ) 1c ( j) yj3hjsn(x ) c (i j) y ,i 0,1,2,L, n（9）ijijn11h(x ) c (n j) ynsnj3j21注意到3 (x) 的局部非零性及其函數(shù)值： 3 (0) 3 ， 3 (1) 6 ，當(dāng) 2 時x 1 ) 知， (0) 0 ， (1) m 1 ， (x) 0 ；且由 (333322 2 時3 (x) 0 。則（9）式的每一個方程只有三個非零系數(shù)，具體為當(dāng) x-187-

38、 c1 c1 2hy0ci1 4ci ci1 6 yi ,i 0,1,2,L, n（10） c c 2hyn1n1n由方程組（10）容易求解出c j（ j 1,0,L, n 1），即表達(dá)式。1.6.4二維等距 B 樣條函數(shù)插值設(shè)有空間曲面 z f (x, y) （未知 ），如 果已知二維等距節(jié)點 (xi , y j ) (x0 ih, y0 j )（ h, 0 ）上的值為 zij （ i 0,1,L, n; j 0,1,L, m ），則相應(yīng)的 B到三次樣條函數(shù) s3 (x)樣條磨光曲面的一般形式為n1 m1 x x0 y y0s(x, y) cijk i l j hi k jl其中 cij （

39、 i 0,1,L, n; j 0,1,L, m ）為待定常數(shù)， k, l 可以取不同值，常用的也是k, l 2 和3 的情形。這是一種具有良好保凸性的光滑曲面（函數(shù)），在工程設(shè)計中是常用的，但只能使用于均勻劃分或近似均勻劃分的情況。1.7 二維插值前面講述的都是一維插值，即節(jié)點為一維變量，插值函數(shù)是一元函數(shù)（曲線）。若節(jié)點是二維的，插值函數(shù)就是二元函數(shù)，即曲面。如在某區(qū)域測量了若干點（節(jié)點）的高程（節(jié)點值），為了畫出較精確的等高線圖，就要先些點的高程（插值）。1.7.1插值節(jié)點為網(wǎng)格節(jié)點的點（插值點），計算這已知 m n 個節(jié)點： (xi , y j , zij ) （ i 1,2,L, m;

40、 j 1,2,L, n ），且 x1 L xm ；y1 L yn 。求點(x, y) 處的插值 z 。中有一些計算二維插值的程序。如z=2(x0,y0,z0,x,y,method)其中 x0,y0 分別為m 維和n 維向量，表示節(jié)點，z0 為n m 維矩陣，表示節(jié)點值，x,y為一維數(shù)組，表示插值點，x 與 y 應(yīng)是方向不同的向量，即一個是行向量，另一個是列向量，z 為矩陣，它的行數(shù)為 y 的維數(shù)，列數(shù)為 x 的維數(shù)，表示得到的插值，method的用法同上面的一維插值。如果是三次樣條插值，可以使用命令pp=cs(x0,y0,z0,conds,valconds)，z=fnval(pp,x,y)其中

41、 x0,y0 分別為m 維和n 維向量，z0 為m n 維矩陣， z 為矩陣，它的行數(shù)為 x 的維數(shù)，列數(shù)為 y 的維數(shù)，表示得到的插值，具體使用方法同一維插值。例2 在一丘陵地帶測量高程，x 和 y 方向每隔100米測一個點，得高程如2表，試插值一曲面，確定合適的模型，并由此找出最高點和該點的高程。解 編寫程序如下：clear,clc x=100:100:500; y=100:100:400;z=636697624478478450420698-188-712630680662pp=cs674626598552412334400310;(x,y,z)xi=100:10:500;yi=100:

42、10:400cz1=fnval(pp,xi,yi)cz2=2(x,y,z,xi,yi,spline)i,j=find(cz1=max(max(cz1) x=xi(i),y=yi(j),zmax=cz1(i,j)表21.7.2插值節(jié)點為散亂節(jié)點已知n 個節(jié)點： (x) ，求點(x, y) 處的插值 z 。對上述問題，中提供了插值函數(shù) griddata，其格式為：ZI = GRIDDATA(X,Y,Z,XI,YI)其中 X、Y、Z 均為n 維向量，指明所給數(shù)據(jù)點的橫坐標(biāo)、縱坐標(biāo)和豎坐標(biāo)。向量 XI、 YI 是給定的網(wǎng)格點的橫坐標(biāo)和縱坐標(biāo)，返回值 ZI 為網(wǎng)格（XI，YI）處的函數(shù)值。XI 與 YI

43、 應(yīng)是方向不同的向量，即一個是行向量，另一個是列向量。例 3在某海域測得一些點（x,y）處的水深 z 由下表給出，在矩形區(qū)域（75,200）(-50,150) 內(nèi)畫出海底曲面的圖形。表 3解 編寫程序如下：x=129 y=7.5z=-4140 103.5141.5 2388 185.5147 22.5195 105 157.5 107.5 7781 162 162 117.5;56.5 -66.5 84 -33.5;137.5 85.5 -6.5 -8138686889988949;xi=75:1:200; yi=-50:1:150;zi=griddata(x,y,z,xi,yi,cubic)

44、subplot(1,2,1), plot(x,y,*)subplot(1,2,2), mesh(xi,yi,zi)2曲線擬合的線性最小二乘法2.1線性最小二乘法曲線擬合問題的提法是，已知一組（二維）數(shù)據(jù)，即平面上的n 個點( xi , yi ) ，i 1,2,L, n ， xi 互不相同，尋求一個函數(shù)（曲線） y f ( x) ，使 f ( x) 在某種準(zhǔn)則下與所有數(shù)據(jù)點最為接近，即曲線擬合得最好。-189-x129140103 588185.5195105157.5107.57781162162117.5y7.5141.52314722.5137.585.5-6.5-81356.5-66.5

45、84-33.5z48686889988949xy100200300400500100636697624478450200698712630478420300680674598412400400662626552334310線性最小二乘法是解決曲線擬合最常用的方法，基本思路是，令f ( x) a1r1 ( x) a2 r2 ( x) L am rm ( x)（11）其中 rk ( x) 是事先選定的一組線性無關(guān)的函數(shù)， ak 是待定系數(shù)(k 1,2,L, m, m n) 。擬合準(zhǔn)則是使 yi ，i 1,2,L, n ，與 f ( xi ) 的距離 i 的平方和最小，稱為最小二乘準(zhǔn)則。2.1.1系

46、數(shù)ak 的確定記nnm iJ (a ,L, a ) 2 f (x ) y 2（12）1iii 1i 1J為求a1 ,L, am 使 J 達(dá)到最小，只需利用極值的必要條件a 0 (k 1,L, m) ，得到k關(guān)于a1,L, am 的線性方程組 rj (xi ) ak rk (xi ) yi 0,( j 1,L, m)nmi1k 1即mnn ak rj (xi )rk (xi ) rj (xi )yi ,( j 1,L, m)（13）k 1i1i1記 r1 ( x1 )Lrm ( x1 )R MMM，r1 ( xn )Lrm ( xn )nmA a1 ,L, am ， Y ( y ,L, y )T

47、T1n方程組（13）可表為RT RA RTY（14）當(dāng)r ( x),L, r ( x) 線性無關(guān)時， R 列滿秩， RT R 可逆，于是方程組（14）有唯1m一解A (RT R)1 RTY2.1.2函數(shù) rk ( x) 的選取面對一組數(shù)據(jù)( x，用線性最小二乘法作曲線擬合時，首要的、也是關(guān)鍵的一步是恰當(dāng)?shù)剡x取 r1( x),L, rm ( x) 。如果通過機(jī)理分析，能夠知道 y 與 x 之樣的函數(shù)關(guān)系，則 r1 ( x),L, rm ( x) 容易確定。若無法知道 y 與 x 之間的關(guān)間應(yīng)該系，通?？梢詫?shù)據(jù)( x擬合。人們常用的曲線有直線y a1x a2多項式y(tǒng) a1x L a x a（一般

48、m 2,3 ，不宜太高）mmm 1作圖，直觀地判斷應(yīng)該用什么樣的曲線去作（iii）雙曲線（一支） y a1 ax-190-2（iv）指數(shù)曲線y a1e 2a x對于指數(shù)曲線，擬合前需作變量代換，化為對a1 , a2 的線性函數(shù)。已知一組數(shù)據(jù)，用什么樣的曲線擬合最好，可以在直觀判斷的基礎(chǔ)上，選幾種曲線分別擬合，然后比較，看哪條曲線的最小二乘指標(biāo) J 最小。最小二乘法的解方程組方法在上面的記號下，實現(xiàn)J (a1,L, am ) | RA Y |2中的線性最小二乘的為RA Y2 ,minA命令為 A R Y 。2例 4用最小二乘法求一個形如 y a bx2 的經(jīng)驗公式，使它與表 4 所示的數(shù)據(jù)擬合。

49、表 4xy192531384419.032.349.073.397.8解 編寫程序如下x=19 y=19.025313844;32.3 49.0 73.397.8;r=ones(5,1),x.2; ab=ryx0=19:0.1:44;y0=ab(1)+ab(2)*x0.2; plot(x,y,o,x0,y0,r)2.2.2多項式擬合方法m ，即用m 次多項式擬合給定數(shù)據(jù)，如果取r ( x),L, r(1m 1中有現(xiàn)成的函數(shù)a=polyfit(x0,y0,m)其中輸入?yún)?shù) x0,y0 為要擬合的數(shù)據(jù)，m 為擬合多項式的次數(shù)，輸出參數(shù) a 為擬合多項式 y=amxm+a1x+a0 系數(shù) a= am

50、,a1, a0。多項式在 x 處的值 y 可用下面的函數(shù)計算 y=polyval(a,x)。例 5 某鄉(xiāng)鎮(zhèn)企業(yè) 1990-1996 年的生產(chǎn)利潤如表 5。表5 年份1990 1991 19921993 1994 1995 1996利潤（萬元） 70122144152174196202試1997 年和 1998 年的利潤。解 作已知數(shù)據(jù)的的散點圖， x0=1990 1991 1992 19931994 1995 1996;y0=70 122 144152 174196 202;plot(x0,y0,*)-191-發(fā)現(xiàn)該鄉(xiāng)鎮(zhèn)企業(yè)的年生產(chǎn)利潤幾乎直線上升。因此，可以用 y a1 x a0 作為擬合函

51、數(shù)來該鄉(xiāng)鎮(zhèn)企業(yè)未來的年利潤。編寫程序如下：x0=1990 1991 1992 1993 1994 1995 1996;y0=70 122 144152 174196 202;a=polyfit(x0,y0,1) y97=polyval(a,1997) y98=polyval(a,1998)求得a1 20 ，a0 4.070510 ，1997 年的生產(chǎn)利潤 y97=233.4286，1998 年的4生產(chǎn)利潤 y98=253.9286。3 最小二乘優(yōu)化在無約束最優(yōu)化問題中，有些重要的特殊情形，比如目標(biāo)函數(shù)由若干個函數(shù)的平方和。這類函數(shù)一般可以寫成：m iF (x) f (x) , x R2ni1n

52、 ) ，一般假設(shè) m n 。Tm其中把極小化這類函數(shù)： iF (x) f (x)2mini1稱為最小二乘優(yōu)化問題。最小二乘優(yōu)化是一類比較特殊的優(yōu)化問題，在處理這類問題時，也提供了優(yōu)化工具箱中，用于求解最小二乘優(yōu)化問題的函數(shù)有：lsqlin、一些強(qiáng)大的函數(shù)。在lsqcurvefit、lsqnonlin、lsqnonneg，用法介紹如下。3.1lsqlin 函數(shù)1222Cx d求解minx A * x bs.t. Aeq * x beqlb x ub其中C, A, Aeq 為矩陣， d , b, beq, lb, ub, x 為向量。中的函數(shù)為：x=lsqlin(C,d,A,b,Aeq,beq,l

53、b,ub,x0)例 6 用lsqlin 命令求解例 4。解 編寫程序如下：x=19 y=19.025313844;32.349.073.397.8;r=ones(5,1),x.2; ab=lsqlin(r,y) x0=19:0.1:44;y0=ab(1)+ab(2)*x0.2; plot(x,y,o,x0,y0,r)-192-3.2 lsqcurvefit 函數(shù)給定輸入輸出數(shù)列 xdata, ydata ,求參量 x ,使得min 1 1 F (x, xdata ) ydata222F (x, xdata) ydataii22xi中的函數(shù)為 X=LSQCURVEFIT(FUN,X0,XDATA

54、,YDATA,LB,UB,OPTIONS)其中 FUN 是定義函數(shù) F (x, xdata) 的 M 文件。例 7 用下面表 6 中的數(shù)據(jù)擬合函數(shù)c(t) a be0 02kt 中的參數(shù) a, b, k 。表6 t j100200300400500600700800900 1000c j4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59解 該問題即解最優(yōu)化問題：min F (a, b, k) (a be100 02ktc )2jji1編寫 M 文件 fun1.m 定義函數(shù) F (x, tdata) ：function f=fun1(x,tdata

55、);f=x(1)+x(2)*exp(-0.02*x(3)*tdata);%其中 x(1)=a，x(2)=b，x(3)=k調(diào)用函數(shù) lsqcurvefit，編寫程序如下： td=100:100:1000;cd=4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59;x0=0.2 0.05 0.05;x=lsqcurvefit(fun1,x0,td,cd)3.3 lsqnonlin 函數(shù)已知函數(shù)向量 F (x) f (x),L,f (x)T ，求 x 使得1kmin 122F (x)2x中的函數(shù)為 X=LSQNONLIN(FUN,X0,LB,UB,OPT

56、IONS)其中 FUN 是定義向量函數(shù) F (x) 的 M 文件。例 8 用lsqnonlin 函數(shù)求解例 7。解 這里F (x) F (x, t) T0 02kta be0 02kt10a bec ,cL,1110 x a,b,k （1）編寫 M 文件 fun2.m 如下： function f=fun2(x); td=100:100:1000;cd=4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59;f=x(1)+x(2)*exp(-0.02*x(3)*td)-cd;-193-（2）調(diào)用函數(shù) lsqnonlin，編寫程序如下： x0=0.2

57、0.05 0.05; %初始值是任意取的 x=lsqnonlin(fun2,x0)3.4 lsqnonneg 函數(shù)1222Cx d求解非負(fù)的 x ，使得滿足 minx中的函數(shù)為，X = LSQNONNE,X0,OPTIONS)0.03720.68610.28690.85870.70710.1781例9 已知C ， d ，求 x(x 0) 滿足0.62330.63440.62450.61700.07470.8405Cx dmin 1222x最小。解 編寫程序如下：c=0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170; d=0.8587

58、;0.1781;0.0747;0.8405;x=lsqnonne)35.曲線擬合的用戶圖形界面求法工具箱提供了命令 cftool，該命令給出了一維數(shù)據(jù)擬合的交互式環(huán)境。具體執(zhí)行步驟如下：把數(shù)據(jù)導(dǎo)入到工作空間；運(yùn)行 cftool，打開用戶圖形界面窗口；對數(shù)據(jù)進(jìn)行預(yù)處理；選擇適當(dāng)?shù)哪Ｐ瓦M(jìn)行擬合；（5）生成一些相關(guān)的統(tǒng)計量，并進(jìn)行?？梢酝ㄟ^幫助（運(yùn)行 doc cftool）熟悉該命令的使用細(xì)節(jié)。4 曲線擬合與函數(shù)近前面講的曲線擬合是已知一組離散數(shù)據(jù)( x，選擇一個較簡單的函數(shù) f ( x) ，如多項式，在一定準(zhǔn)則如最小二乘準(zhǔn)則下，最接近這些數(shù)據(jù)。如果已知一個較為復(fù)雜的連續(xù)函數(shù) y( x), x a,

59、 b，要求選擇一個較簡單的函數(shù)f ( x) ，在一定準(zhǔn)則下最接近 f ( x) ，就是所謂函數(shù)近。與曲線擬合的最小二乘準(zhǔn)則相對應(yīng)，函數(shù)近常用的一種準(zhǔn)則是最小平方近，即bJ f ( x) y( x)2 dx(15)a達(dá)到最小。與曲線擬合一樣，選一組函數(shù)rk ( x), k 1,L, m構(gòu)造 f ( x) ，即令f ( x) a1r1 ( x) L am rm ( x)-194-代入（15）式，求a1 ,L, am 使 J 達(dá)到極小。利用極值必要條件 (r1, r1 )(r1, rm ) a1 ( y, r1 ) L M MMMM（16） (rm , r1 )L(rm , rm ) am ( y,

60、 rm )b這里(g, h) g( x)h( x)dx 。當(dāng)方程組（16）的系數(shù)矩陣非奇異時，有唯一解。a最簡單的當(dāng)然是用多項式近函數(shù)，即選 r ( x) 1 ， r ( x) x ， r ( x) x2 ，。123b( x)dx 0 ，(i j) ，方程組（16）的系數(shù)矩陣將是對角陣，計并且如果能使 r ( x)rija算大大簡化。滿足這種性質(zhì)的多項式稱正交多項式。勒讓得（Legendre）多項式是在1,1 區(qū)間上的正交多項式，它的表達(dá)式為d k1P0 ( x) 1,Pk ( x) 2k k! dxk ( x 1) , k 1,2,L2k可以證明i j0,P ( x)P ( x)dx 12i