1、第四講 插值與擬合之插值(上),內容：插值是離散函數逼近的重要方法，利用它 可通過函數在有限個點處的取值狀況，估 算出函數在其他點處的近似值 目的：學習插值的基本思想和方法，掌握Matlab 的一維/二維等距和非等距插值函數 要求：掌握Matlab插值函數，處理插值應用問題 了解拉格朗日和分段線性插值的基本思想 了解三次樣條插值的提法和思路 掌握插值函數 interp interp1 interp2 griddata 掌握水塔用水量的計算(水位-體積-流速-積分),關于插值與擬合的區(qū)別,面對工程實踐和科學計算中的采集得到數據(xi,yi)，我們總是試圖去揭示x與y之間的關系，即用近似的y=f(

2、x)來表示，那么我們通常可以采用兩種方法：插值與擬合 插值與擬合的區(qū)別在于 插值試圖去通過已知點了解未知 點處的函數值；而擬合則在于在 整體上用某種已知函數去擬合數 據點列所在未知函數的性態(tài)。 關鍵區(qū)別在于插值要求必須經過已知點列，擬合只求盡量靠近不必經過！擬合將在本講下介紹,引例1 函數查表問題： 已知標準正態(tài)分布函數表，求表中沒有的值 (2.34)=0.99036 (2.35)=0.99061 求 (2.3457) (2.35-2.3457)/(2.35-2.34)* (2.34)+ (2.3457-2.34)/(2.35-2.34)* (2.35) 引例2 地圖繪制問題： 假如我們在地圖

3、邊界獲取了一些邊界點的坐標，連接這些邊界點形成閉合曲線，可以用來近似表示真實邊界線，如何更準確地逼近真實邊界線？,函數查表與地圖邊界線繪制,如何更準確地逼近真實邊界線？,插值在數碼圖像放大中的應用,引例3 圖像插值放大： 數碼相機運用插值的方法可以創(chuàng)造出比傳感器實際像素更多的圖像，這種處理稱為“數碼變焦”。,106*40原始圖像：,左邊： 最近鄰插值 放大450%,右邊： 雙三次插值 放大450%,插值在圖像三維重建中的應用,Surface recostruction from scattered points cloud,分段線性插值和拉格朗日插值,分段線性插值： 用直線(線性)連接數據點列

4、上相鄰的兩點。 比如在兩點xi-1,xi上線性插值函數為,拉格朗日插值： 用n次拉格朗日插值多項式,連接數據點列上相鄰的n+1個點。Pszjs71,拉格朗日插值基函數的構造,比如 在三個點x0,x1,x2上lagrange插值函數為 (線性插值是拉格朗日插值最簡單的情形),分段三次埃爾米特插值條件數,分段三次埃爾米特插值： 線性插值在每一小段上(兩點之間)，用到2個條件q(xi)=yi，所以確定了一個線性插值函數；三次埃爾米特插值在每一小段上，用到4個條件q(xi)=yi， q(xi)=yi,所以確定一個3次多項式插值函數。 分段插值主要是為了避免高次插值可能出現(xiàn)的大幅度振蕩現(xiàn)象，在實際應用中

5、通常采用分段低次插值來提高近似程度，比如可用分段線性插值或分段三次埃爾米特插值來逼近已知函數，但它們的總體光滑性較差，為了克服這一缺點，三次樣條插值成為比較理想的工具。,三次樣條(spline)插值的概念,樣條的概念出自工程設計和機械加工(飛機、船舶外形曲線設計)中的繪圖工具(曲線尺)，簡單說就是具有連續(xù)二階導數的三次插值多項式函數。,三次樣條(spline)插值的條件數,首先從段數n=2分析：我們知道在每一小段的三次多項式有4個系數，所以如下圖，總共需要有4*2=8個方程來確定； 由q(xi)=yi可以確定2*2=4個方程，又由內部節(jié)點q1(xi)= q2(xi)和q1(xi)= q2(xi

6、)可以確定2*(2-1)=2個方程，看來剩下的8-(4+2)=2個方程只有靠外部給定(邊界條件)了,一維曲線等距插值函數interp,interps syntax One-dimensional r times longer data interpolation y = interp(y,r) 題例 在原始數據點中增倍插值 x=0:0.001:1; y=sin(2*pi*30*x)+sin(2*pi*60*x); yi=interp(y,4); subplot(1,2,1); stem(y(1:30); title(Original Points); subplot(1,2,2); stem(

7、yi(1:120); title(Interpolated Points);,一維曲線等距插值函數interp1,interp1s syntax One-dimensional data interpolation yi = interp1(x,y,xi,method) nearest Nearest neighbor interpolation linear Linear interpolation (default) spline Cubic spline interpolation cubic Piecewise cubic Hermite interpolation 題例 在一天24小

8、時內,從零點開始每間隔2小時測得的環(huán)境溫度，推測在15點6分的的溫度 x=0:2:24; y=12,9,9,10,18,24,28,27,25,20,18,15,13; plot(x,y,-ro); hold on; xi=15.1; yi=interp1(x,y,xi,spline), xi=0:1/3600:24; yi=interp1(x,y,xi,spline); plot(xi,yi,b-);,二維曲面等距插值函數interp2,interp2s syntax Two-dimensional data interpolation ZI = interp2(X,Y,Z,XI,YI,me

9、thod) nearest Nearest neighbor interpolation linear Bilinear interpolation (default) spline Cubic spline interpolation cubic Bicubuc interpolation,二維曲面等距插值函數interp2,動畫展示：三維空間中的曲面等距格點,二維曲面等距插值函數interp2,題例 粗糙山頂曲面的平滑處理(等距情形) load mountain.mat %載入山頂地形數據 mesh(x,y,z) %繪制原始山頂地形圖 xi=linspace(0,5,50); yi=lin

10、space(0,6,80); xii,yii=meshgrid(xi,yi); zii=interp2(x,y,z,xii,yii,spline); %三次樣條插值 figure; surf(xii,yii,zii) %繪制平滑處理后的山頂曲面 hold on; xx,yy=meshgrid(x,y); plot3(xx,yy,z+0.1,ob);,二維曲面等距插值函數interp2,題例 粗糙山頂曲面的平滑處理(等距情形),二維曲面散亂插值函數griddata,griddatas syntax Data interpolation for scattered points ZI = grid

11、data(x,y,z,XI,YI) XI,YI,ZI = griddata(x,y,z,xi,yi) . = griddata(.,method) linear Triangle-based linear interpolation cubic Triangle-based cubic (default) nearest Nearest neighbor v4 MATLAB 4 griddata method MATLAB二維插值函數griddata， 可以將平面或曲面上的散亂點插值為規(guī)則網格,二維曲面散亂插值函數griddata,題例 粗糙山頂曲面的平滑處理(散亂情形) rand(seed,

12、0) x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(-x.2-y.2); plot3(x,y,z,o); hold on ti = -2:.25:2; XI,YI = meshgrid(ti,ti); ZI = griddata(x,y,z,XI,YI); mesh(XI,YI,ZI);,二維曲面散亂插值函數griddata,題例 墨西哥草帽的平滑處理(散亂情形) x = rand(100,1)*16 - 8; y = rand(100,1)*16 - 8; r = sqrt(x.2 + y.2) + eps; z = sin(r).

13、/r; plot3(x,y,z,.,MarkerSize,15) hold on xlin = linspace(min(x),max(x),33); ylin = linspace(min(y),max(y),33); X,Y = meshgrid(xlin,ylin); Z = griddata(x,y,z,X,Y,cubic); mesh(X,Y,Z); axis tight;,南半球氣旋變化的可視圖形,山區(qū)地貌的可視化圖形,水塔用水量估計通用程序,通用程序tbp69.m可近似計算時間段內的用水量 格式為：tbp69(ts,tf) 其中ts為起點時間，tf為終點時間,實驗一：水塔用水量估

14、計Thats all3Q!,第四講 插值與擬合之擬合(下),內容：擬合是離散函數逼近的重要方法，利用它 可通過函數在有限個點處的取值狀況，擬 合出近似替代函數，進而估算出函數在其 他點處的近似值。 目的：學習擬合的基本思想和方法，掌握Matlab 的多項式/一般擬合函數/曲線擬合工具箱 要求：掌握Matlab擬合函數，處理擬合應用問題 了解基于最小二乘法則擬合的基本思想 掌握擬合函數 polyfit lsqcurvefit curvefit 掌握cftool曲線擬合工具箱(多目標函數多法則),關于數據擬合的兩個要素.,在工程實踐和科學計算中，用某種經驗函數解析式y(tǒng)=f(x)來近似刻畫采集數據(

15、x,y) 之間的關系的方法就叫擬合，所謂“擬合”有 “最貼近”之意 。 與插值不同，擬合的主要目 標是要離散點盡量靠近擬合函 數。一般過程是，我們首先根 據采樣點的散點分布圖，大致 推測x與y之間的經驗函數形式 (比如多項式、指數函數等)， 然后依據某種法則(比如最常用的最小二乘法則)，確定出的經驗函數解析式中的待定參數。其中經驗函數和擬合法則是擬合的兩個關鍵要素！,引例 1 化合物濃度隨時間變化的規(guī)律： 與插值面臨的問題相似，我們被要求去求解或預測表格中沒有的因變量對應值，與插值的解決思路不同，我們試圖獲得比較完備的解決方案：設計并求出離散數據點的近似替代函數，有了近似函數解析式，就可以進一

16、步代值計算或作圖分析。 為了揭示濃度y與時間t之間呈現(xiàn)的函數規(guī)律，我們首先作出散點圖，幫助分析和設計經驗函數,化合物濃度隨時間變化的規(guī)律,化合物濃度隨時間變化的規(guī)律,如圖, 化合物濃度y隨時間t大致呈拋物線狀（二次函數）變化，這種分析和判斷來自已有經驗.,t=1:16; c=4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6; plot(t,c,-ro),化合物濃度隨時間變化的規(guī)律,經驗函數形式： 已經擬定為多項式函數：y= at2 +bt+ c 剩下的工作是確定擬合原則： 可選的法則很多，其中最常用的

17、是最小二乘法則(method of Least Squares)，即各點偏差平方和最小,高斯和勒讓德關于最小二乘法的發(fā)明權,化合物濃度隨時間變化的規(guī)律,對經驗函數形式確定的補充說明： 擬合函數解析式選用什么形式(用多項式，還是用冪函數?)，主要取決于采樣點的分布無疑，那么如何求出這些含有待定參數的解析式呢？把各點偏差的平方和最小作為一個目標函數，實際上考慮為極值問題，極值點導數為零。具體計算時，我們在把經驗函數用一系列擬合基函數線性表出同時，在j個采樣點對待定參數Cj求偏導(=0),獲取j個方程，進而解出Cj，具體參見數值計算SZJSp9091 在本例中，已經擬定擬合的目標函數為多項式函數：y

18、= at2 +bt+ c ，所以只要解出三個待定參數a,b,c，問題即獲解決,基于最小二乘的多項式擬合函數,基于最小二乘的多項式擬合函數polyfit： Polynomial curve fitting .Syntax： p = polyfit(x,y,n) 其中n是擬合多項式的階數，不能超過（散點數據對數-1） 下面回到化合物濃度隨時間變化的引例： t=1:16; c=4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6; plot(t,c,ko); hold on; %作散點圖 p2=polyfit(

19、t,c,2); y2=poly2str(p2,t), %作多次擬合比較 p5=polyfit(t,c,5); y5=poly2sym(p5,t), f=inline(y5) ti=0:.001:20; plot(ti,polyval(p2,ti),b-,ti,f(ti),r-); disp(化合物在刻度11.2的濃度近似值為,num2str(f(11.2) disp(化合物在刻度17.8的濃度預測值為,num2str(f(17.8) stem(11.2 17.8,f(11.2) f(17.8),r); xlabel(時間t); ylabel(化合物濃度c); title(化合物濃度隨時間變化的

20、規(guī)律),引例 2 確定醫(yī)用薄膜滲透率的數學模型： 某種醫(yī)用薄膜允許一種物質分子 從高濃度溶液VB穿過薄膜向低濃度 溶液VA中擴散。通過單位面積膜S 分子擴散的速度與膜兩側溶液的濃 度差成正比，比例系數K表示薄膜 被該物質分子穿透的能力，稱為滲透率，定時測量薄膜VB側的溶液濃度值CB，以此確定K的值VA=VB=1000cm3, S=10cm2, 容器的B部分溶液濃度CB的測試結果如下表：( CB單位為mg/ cm3 ),確定醫(yī)用薄膜滲透率的數學模型,確定醫(yī)用薄膜滲透率的數學模型,由質量守恒考察 t,t+t 時間段B向A中滲透物質： VA*CA(t+t)-VA*CA(t) = SKCB(t)-CA

21、(t)t 推出 dCA(t)/dt = SK/VA*CB(t)-CA(t) 兩邊除以t, 令t0 又由質量守恒考察整個容器中物質總量始終不變： VA*CA(t)+VB*CB(t) = VA*aA+VB*aB 推出 CA(t) = aA+VB/VA*aB-VB/VA*CB(t) 代入上式2 推出 dCB(t)/dt = -SK(1/VA+1/VB)CB(t)+SK(aA/VB+aB/VA) CB(0)=aB 初值條件 此帶初值微分方程可由dsolve求解,在上式中，已知的包括VA，VB，S以及一組t和CB(t)值 未知的包括aA，aB，K，下面通過數據擬合確定滲透率K,確定醫(yī)用薄膜滲透率的數學模

22、型,在上式中，代入已知值VA=VB=1000cm3，S=10cm2 令a=(aA*VA+aB*VB)/(VA+VB)，b=VA(aB-aA)/(VA+VB) 簡化之后的表達式為：CB(t)=a+b*exp(-0.02*k*t) 編寫被調M文件 tbp79.m function CB=tbp79(x,t) CB=x(1)+x(2)*exp(-0.02*x(3)*t); 編寫主調M文件 fittbp79.m（片段） x=curvefit(tbp79,x0,t,CB) %curvefit擬合及圖像 x=lsqcurvefit(tbp79,x0,t,CB) %lsqcurvefit擬合及圖像 求解結果

23、：a=x(1)=0.0070; b=x(2)=-0.0030; k=x(3)=0.1012 進一步求解：aA=0.01；aB=0.004 最終數學模型：CB(t)=0.007-0.003*exp(-0.002*t),基于最小二乘的一般擬合函數,基于最小二乘的一般擬合函數lsqcurvefit： Solve nonlinear curve-fitting (data-fitting) problems in the least-squares sense.Syntax： x = lsqcurvefit(fun,x0,xdata,ydata) x,resnorm = lsqcurvefit(.)

24、范例: function F=myfun(x,xdata) F=x(1)*xdata.2+x(2)*sin(xdata)+x(3)*xdata.3; %下面是主調函數fitmyfun.m的部分 xdata = 3.6 7.7 9.3 4.1 8.6 2.8 1.3 7.9 10.0 5.4; ydata = 16.5 150.6 263.1 24.7 208.5 9.9 2.7 163.9 325.0 54.3; x0 = 10, 10, 10;% Starting guess x,resnorm = lsqcurvefit(myfun,x0,xdata,ydata),MATLAB工具箱的版本

25、更新.,CFTool曲線擬合工具箱簡介,基于MATLAB的曲線擬合問題，已經提供獨立的toolbox供調用，該toolbox采用GUI界面，功能強大，下面簡單介紹如何使用該Toolbox解決一般曲線擬合問題。 在command window中鍵入指令cftool即可啟動曲線擬合工具箱。在該集成環(huán)境里面，可以實現(xiàn)多種經驗函數，多種法則的曲線擬合，實時繪制圖像并進行誤差分析。 需要注意的是：在進入Curve Fitting Toolbox環(huán)境進行曲線擬合之前，需要預先在workspace輸入或載入供擬合的數據源,CFTool-選擇Data導入數據,下面還是以引例的采樣數據為例，進行演示： t=1:

26、16;y=4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6;cftool,導入數據,繪制散點圖,CFTool-選擇Fitting擬合數據,進行擬合,這里可供選擇的擬合類型和可選參數比較多,包括多項式函數,指數函數,冪函數等,如何確定最優(yōu)的方案？,CFTool-擬合效果評價指標,SSE - The sum of squares due to error. This statistic measures the deviation of the responses from the fitted values of the responses.