曲線曲面的插值與擬合方法.ppt

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

2、x)來表示，那么我們通?？梢圆捎脙煞N方法：插值與擬合 插值與擬合的區(qū)別在于 插值試圖去通過已知點(diǎn)了解未知 點(diǎn)處的函數(shù)值；而擬合則在于在 整體上用某種已知函數(shù)去擬合數(shù) 據(jù)點(diǎn)列所在未知函數(shù)的性態(tài)。 關(guān)鍵區(qū)別在于插值要求必須經(jīng)過已知點(diǎn)列，擬合只求盡量靠近不必經(jīng)過！擬合將在本講下介紹,引例1 函數(shù)查表問題： 已知標(biāo)準(zhǔn)正態(tài)分布函數(shù)表，求表中沒有的值 (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、邊界獲取了一些邊界點(diǎn)的坐標(biāo)，連接這些邊界點(diǎn)形成閉合曲線，可以用來近似表示真實(shí)邊界線，如何更準(zhǔn)確地逼近真實(shí)邊界線？,函數(shù)查表與地圖邊界線繪制,如何更準(zhǔn)確地逼近真實(shí)邊界線？,插值在數(shù)碼圖像放大中的應(yīng)用,引例3 圖像插值放大： 數(shù)碼相機(jī)運(yùn)用插值的方法可以創(chuàng)造出比傳感器實(shí)際像素更多的圖像，這種處理稱為“數(shù)碼變焦”。,106*40原始圖像：,左邊： 最近鄰插值 放大450%,右邊： 雙三次插值 放大450%,插值在圖像三維重建中的應(yīng)用,Surface recostruction from scattered points cloud,分段線性插值和拉格朗日插值,分段線性插值： 用直線(線性)連接數(shù)據(jù)點(diǎn)列

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

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

6、)可以確定2*(2-1)=2個方程，看來剩下的8-(4+2)=2個方程只有靠外部給定(邊界條件)了,一維曲線等距插值函數(shù)interp,interps syntax One-dimensional r times longer data interpolation y = interp(y,r) 題例 在原始數(shù)據(jù)點(diǎn)中增倍插值 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);,一維曲線等距插值函數(shù)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、時內(nèi),從零點(diǎn)開始每間隔2小時測得的環(huán)境溫度，推測在15點(diǎn)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-);,二維曲面等距插值函數(shù)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,二維曲面等距插值函數(shù)interp2,動畫展示：三維空間中的曲面等距格點(diǎn),二維曲面等距插值函數(shù)interp2,題例 粗糙山頂曲面的平滑處理(等距情形) load mountain.mat %載入山頂?shù)匦螖?shù)據(jù) mesh(x,y,z) %繪制原始山頂?shù)匦螆D 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);,二維曲面等距插值函數(shù)interp2,題例 粗糙山頂曲面的平滑處理(等距情形),二維曲面散亂插值函數(shù)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二維插值函數(shù)griddata， 可以將平面或曲面上的散亂點(diǎn)插值為規(guī)則網(wǎng)格,二維曲面散亂插值函數(shù)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);,二維曲面散亂插值函數(shù)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ū)地貌的可視化圖形,水塔用水量估計(jì)通用程序,通用程序tbp69.m可近似計(jì)算時間段內(nèi)的用水量 格式為：tbp69(ts,tf) 其中ts為起點(diǎn)時間，tf為終點(diǎn)時間,實(shí)驗(yàn)一：水塔用水量估

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

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

16、步代值計(jì)算或作圖分析。 為了揭示濃度y與時間t之間呈現(xiàn)的函數(shù)規(guī)律，我們首先作出散點(diǎn)圖，幫助分析和設(shè)計(jì)經(jīng)驗(yàn)函數(shù),化合物濃度隨時間變化的規(guī)律,化合物濃度隨時間變化的規(guī)律,如圖, 化合物濃度y隨時間t大致呈拋物線狀（二次函數(shù)）變化，這種分析和判斷來自已有經(jīng)驗(yàn).,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ī)律,經(jīng)驗(yàn)函數(shù)形式： 已經(jīng)擬定為多項(xiàng)式函數(shù)：y= at2 +bt+ c 剩下的工作是確定擬合原則： 可選的法則很多，其中最常用的

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

18、= at2 +bt+ c ，所以只要解出三個待定參數(shù)a,b,c，問題即獲解決,基于最小二乘的多項(xiàng)式擬合函數(shù),基于最小二乘的多項(xiàng)式擬合函數(shù)polyfit： Polynomial curve fitting .Syntax： p = polyfit(x,y,n) 其中n是擬合多項(xiàng)式的階數(shù)，不能超過（散點(diǎn)數(shù)據(jù)對數(shù)-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; %作散點(diǎn)圖 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的濃度預(yù)測值為,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ī)用薄膜滲透率的數(shù)學(xué)模型： 某種醫(yī)用薄膜允許一種物質(zhì)分子 從高濃度溶液VB穿過薄膜向低濃度 溶液VA中擴(kuò)散。通過單位面積膜S 分子擴(kuò)散的速度與膜兩側(cè)溶液的濃 度差成正比，比例系數(shù)K表示薄膜 被該物質(zhì)分子穿透的能力，稱為滲透率，定時測量薄膜VB側(cè)的溶液濃度值CB，以此確定K的值VA=VB=1000cm3, S=10cm2, 容器的B部分溶液濃度CB的測試結(jié)果如下表：( CB單位為mg/ cm3 ),確定醫(yī)用薄膜滲透率的數(shù)學(xué)模型,確定醫(yī)用薄膜滲透率的數(shù)學(xué)模型,由質(zhì)量守恒考察 t,t+t 時間段B向A中滲透物質(zhì)： VA*CA(t+t)-VA*CA(t) = SKCB(t)-CA

21、(t)t 推出 dCA(t)/dt = SK/VA*CB(t)-CA(t) 兩邊除以t, 令t0 又由質(zhì)量守恒考察整個容器中物質(zhì)總量始終不變： 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，下面通過數(shù)據(jù)擬合確定滲透率K,確定醫(yī)用薄膜滲透率的數(shù)學(xué)模

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

23、：a=x(1)=0.0070; b=x(2)=-0.0030; k=x(3)=0.1012 進(jìn)一步求解：aA=0.01；aB=0.004 最終數(shù)學(xué)模型：CB(t)=0.007-0.003*exp(-0.002*t),基于最小二乘的一般擬合函數(shù),基于最小二乘的一般擬合函數(shù)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; %下面是主調(diào)函數(shù)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的曲線擬合問題，已經(jīng)提供獨(dú)立的toolbox供調(diào)用，該toolbox采用GUI界面，功能強(qiáng)大，下面簡單介紹如何使用該Toolbox解決一般曲線擬合問題。 在command window中鍵入指令cftool即可啟動曲線擬合工具箱。在該集成環(huán)境里面，可以實(shí)現(xiàn)多種經(jīng)驗(yàn)函數(shù)，多種法則的曲線擬合，實(shí)時繪制圖像并進(jìn)行誤差分析。 需要注意的是：在進(jìn)入Curve Fitting Toolbox環(huán)境進(jìn)行曲線擬合之前，需要預(yù)先在workspace輸入或載入供擬合的數(shù)據(jù)源,CFTool-選擇Data導(dǎo)入數(shù)據(jù),下面還是以引例的采樣數(shù)據(jù)為例，進(jìn)行演示： 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,導(dǎo)入數(shù)據(jù),繪制散點(diǎn)圖,CFTool-選擇Fitting擬合數(shù)據(jù),進(jìn)行擬合,這里可供選擇的擬合類型和可選參數(shù)比較多,包括多項(xiàng)式函數(shù),指數(shù)函數(shù),冪函數(shù)等,如何確定最優(yōu)的方案？,CFTool-擬合效果評價指標(biāo),SSE - The sum of squares due to error. This statistic measures the deviation of the responses from the fitted values of the responses.