四講曲線曲面插值與擬合方法次課

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

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

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

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

6、有靠外部給定(邊界條件)了q1q2x0 x1x2一維曲線等距插值函數(shù)interpinterps syntaxOne-dimensional r times longer data interpolation y = interp(y,r)題例 在原始數(shù)據(jù)點中增倍插值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(yi(1:120); title(Interpolated Po

7、ints);一維曲線等距插值函數(shù)interp1interp1s syntaxOne-dimensional data interpolation yi = interp1(x,y,xi,method)nearest Nearest neighbor interpolationlinear Linear interpolation (default)spline Cubic spline interpolationcubic Piecewise cubic Hermite interpolation題例 在一天24小時內,從零點開始每間隔2小時測得的環(huán)境溫度，推測在15點6分的的溫度x=0:2:

8、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ù)interp2interp2s syntaxTwo-dimensional data interpolation ZI = interp2(X,Y,Z,XI,YI,method)nearest Nearest neighbor interpolationlin

9、ear Bilinear interpolation (default)spline Cubic spline interpolationcubic Bicubuc interpolation二維曲面等距插值函數(shù)interp2動畫展示：三維空間中的曲面等距格點二維曲面等距插值函數(shù)interp2題例 粗糙山頂曲面的平滑處理(等距情形)load mountain.mat %載入山頂?shù)匦螖?shù)據(jù)mesh(x,y,z) %繪制原始山頂?shù)匦螆Dxi=linspace(0,5,50); yi=linspace(0,6,80); xii,yii=meshgrid(xi,yi); zii=interp2(x,y,z

10、,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ù)griddatagriddatas syntaxData interpolation for scattered points ZI = griddata(x,y,z,XI,YI)XI,YI,ZI = griddata(x,y,z,xi,yi). = griddata(.,met

11、hod)linear Triangle-based linear interpolationcubic Triangle-based cubic (default) nearest Nearest neighbor v4 MATLAB 4 griddata methodMATLAB二維插值函數(shù)griddata，可以將平面或曲面上的散亂點插值為規(guī)則網(wǎng)格二維曲面散亂插值函數(shù)griddata題例 粗糙山頂曲面的平滑處理(散亂情形)rand(seed,0)x = rand(100,1)*4-2; y = rand(100,1)*4-2;z = x.*exp(-x.2-y.2);plot3(x,y,z,

12、o);hold onti = -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)./r;plot3(x,y,z,.,MarkerSize,15)hold onxlin = linspace(min(x),max(x),33);ylin = linspace(

13、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為終點時間實驗一：水塔用水量估計Thats all3Q!第四講 插值與擬合之擬合(下)內容：擬合是離散函數(shù)逼近的重要方法，利用它 可通過函數(shù)在有限個點處的取值狀況，擬 合出近似替代函數(shù)，進而估算出函數(shù)在其 他點處的近似值。目的：學

14、習擬合的基本思想和方法，掌握Matlab 的多項式/一般擬合函數(shù)/曲線擬合工具箱要求：掌握Matlab擬合函數(shù)，處理擬合應用問題了解基于最小二乘法則擬合的基本思想掌握擬合函數(shù) polyfit lsqcurvefit curvefit掌握cftool曲線擬合工具箱(多目標函數(shù)多法則)關于數(shù)據(jù)擬合的兩個要素. 在工程實踐和科學計算中，用某種經(jīng)驗函數(shù)解析式y(tǒng)=f(x)來近似刻畫采集數(shù)據(jù)(x,y) 之間的關系的方法就叫擬合，所謂“擬合”有 “最貼近”之意 。 與插值不同，擬合的主要目標是要離散點盡量靠近擬合函數(shù)。一般過程是，我們首先根據(jù)采樣點的散點分布圖，大致推測x與y之間的經(jīng)驗函數(shù)形式(比如多項式、

15、指數(shù)函數(shù)等)， 然后依據(jù)某種法則(比如最常用的最小二乘法則)，確定出的經(jīng)驗函數(shù)解析式中的待定參數(shù)。其中經(jīng)驗函數(shù)和擬合法則是擬合的兩個關鍵要素！ 引例 1 化合物濃度隨時間變化的規(guī)律： 與插值面臨的問題相似，我們被要求去求解或預測表格中沒有的因變量對應值，與插值的解決思路不同，我們試圖獲得比較完備的解決方案：設計并求出離散數(shù)據(jù)點的近似替代函數(shù)，有了近似函數(shù)解析式，就可以進一步代值計算或作圖分析。 為了揭示濃度y與時間t之間呈現(xiàn)的函數(shù)規(guī)律，我們首先作出散點圖，幫助分析和設計經(jīng)驗函數(shù)化合物濃度隨時間變化的規(guī)律化合物濃度隨時間變化的規(guī)律如圖, 化合物濃度y隨時間t大致呈拋物線狀（二次函數(shù)）變化，這種分

16、析和判斷來自已有經(jīng)驗.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)驗函數(shù)形式：已經(jīng)擬定為多項式函數(shù)：y= at2 +bt+ c剩下的工作是確定擬合原則：可選的法則很多，其中最常用的是最小二乘法則(method of Least Squares)，即各點偏差平方和最小高斯和勒讓德關于最小二乘法的發(fā)明權化合物濃度隨時間變化的規(guī)律對經(jīng)驗函數(shù)形式確定的補充說明： 擬合函數(shù)解析式選用什么形式(用多項式，還是用冪函數(shù)?)，主要取決于采樣

17、點的分布無疑，那么如何求出這些含有待定參數(shù)的解析式呢？把各點偏差的平方和最小作為一個目標函數(shù)，實際上考慮為極值問題，極值點導數(shù)為零。具體計算時，我們在把經(jīng)驗函數(shù)用一系列擬合基函數(shù)線性表出同時，在j個采樣點對待定參數(shù)Cj求偏導(=0),獲取j個方程，進而解出Cj，具體參見數(shù)值計算SZJSp9091 在本例中，已經(jīng)擬定擬合的目標函數(shù)為多項式函數(shù)：y= at2 +bt+ c ，所以只要解出三個待定參數(shù)a,b,c，問題即獲解決基于最小二乘的多項式擬合函數(shù)基于最小二乘的多項式擬合函數(shù)polyfit：Polynomial curve fitting .Syntax：p = polyfit(x,y,n)其中

18、n是擬合多項式的階數(shù)，不能超過（散點數(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; %作散點圖p2=polyfit(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

19、-); 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(化合物濃度隨時間變化的規(guī)律)引例 2 確定醫(yī)用薄膜滲透率的數(shù)學模型： 某種醫(yī)用薄膜允許一種物質分子從高濃度溶液VB穿過薄膜向低濃度溶液VA中擴散。通過單位面積膜S分子擴散的速度與膜兩側溶液的濃度差成正比，比例系數(shù)K表示薄膜被該物質分子穿透的能力，稱為滲透率，定時測量薄膜VB側的溶液濃度值CB，以此確

20、定K的值VA=VB=1000cm3, S=10cm2, 容器的B部分溶液濃度CB的測試結果如下表：( CB單位為mg/ cm3 )確定醫(yī)用薄膜滲透率的數(shù)學模型確定醫(yī)用薄膜滲透率的數(shù)學模型由質量守恒考察 t,t+t 時間段B向A中滲透物質：VA*CA(t+t)-VA*CA(t) = SKCB(t)-CA(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)

21、/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ù)學模型在上式中，代入已知值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(

22、2)*exp(-0.02*x(3)*t);編寫主調M文件 fittbp79.m（片段）x=curvefit(tbp79,x0,t,CB) %curvefit擬合及圖像x=lsqcurvefit(tbp79,x0,t,CB) %lsqcurvefit擬合及圖像求解結果：a=x(1)=0.0070; b=x(2)=-0.0030; k=x(3)=0.1012進一步求解：aA=0.01；aB=0.004最終數(shù)學模型：CB(t)=0.007-0.003*exp(-0.002*t)基于最小二乘的一般擬合函數(shù)基于最小二乘的一般擬合函數(shù)lsqcurvefit：Solve nonlinear curve-fi

23、tting (data-fitting) problems in the least-squares sense.Syntax：x = lsqcurvefit(fun,x0,xdata,ydata)x,resnorm = lsqcurvefit(.)范例:function F=myfun(x,xdata)F=x(1)*xdata.2+x(2)*sin(xdata)+x(3)*xdata.3;%下面是主調函數(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

24、.5 9.9 2.7 163.9 325.0 54.3;x0 = 10, 10, 10;% Starting guessx,resnorm = lsqcurvefit(myfun,x0,xdata,ydata)MATLAB工具箱的版本更新. 同名升級換名升級CFTool曲線擬合工具箱簡介 基于MATLAB的曲線擬合問題，已經(jīng)提供獨立的toolbox供調用，該toolbox采用GUI界面，功能強大，下面簡單介紹如何使用該Toolbox解決一般曲線擬合問題。 在command window中鍵入指令cftool即可啟動曲線擬合工具箱。在該集成環(huán)境里面，可以實現(xiàn)多種經(jīng)驗函數(shù)，多種法則的曲線擬合，實時

25、繪制圖像并進行誤差分析。 需要注意的是：在進入Curve Fitting Toolbox環(huán)境進行曲線擬合之前，需要預先在workspace輸入或載入供擬合的數(shù)據(jù)源CFTool-選擇Data導入數(shù)據(jù)下面還是以引例的采樣數(shù)據(jù)為例，進行演示：t=1: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導入數(shù)據(jù)繪制散點圖CFTool-選擇Fitting擬合數(shù)據(jù)進行擬合擬合方法結果和誤差分析 這里可供選擇的擬合類型和可選參數(shù)比較多,包括多項式函數(shù),指數(shù)函數(shù),冪函數(shù)等,如何確定最優(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. A