第四講曲線曲面的插值與擬合方法次課演示文稿

第四講曲線曲面的插值與擬合方法次課演示文稿當(dāng)前1頁，總共42頁。優(yōu)選第四講曲線曲面的插值與擬合方法次課當(dāng)前2頁，總共42頁。引例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地圖繪制問題：假如我們在地圖邊界獲取了一些邊界點(diǎn)的坐標(biāo)，連接這些邊界點(diǎn)形成閉合曲線，可以用來近似表示真實(shí)邊界線，如何更準(zhǔn)確地逼近真實(shí)邊界線？函數(shù)查表與地圖邊界線繪制(2.3457)=？當(dāng)前3頁，總共42頁。如何更準(zhǔn)確地逼近真實(shí)邊界線？當(dāng)前4頁，總共42頁。插值在數(shù)碼圖像放大中的應(yīng)用引例3圖像插值放大：數(shù)碼相機(jī)運(yùn)用插值的方法可以創(chuàng)造出比傳感器實(shí)際像素更多的圖像，這種處理稱為“數(shù)碼變焦”。106*40原始圖像：左邊：最近鄰插值放大450%右邊：雙三次插值放大450%當(dāng)前5頁，總共42頁。插值在圖像三維重建中的應(yīng)用Surfacerecostructionfromscatteredpointscloud

當(dāng)前6頁，總共42頁。分段線性插值和拉格朗日插值分段線性插值：用直線(線性)連接數(shù)據(jù)點(diǎn)列上相鄰的兩點(diǎn)。比如~在兩點(diǎn)[xi-1,xi]上線性插值函數(shù)為~拉格朗日插值：用n次拉格朗日插值多項(xiàng)式連接數(shù)據(jù)點(diǎn)列上相鄰的n+1個(gè)點(diǎn)。Pszjs71當(dāng)前7頁，總共42頁。拉格朗日插值基函數(shù)的構(gòu)造比如在三個(gè)點(diǎn)[x0,x1,x2]上lagrange插值函數(shù)為(線性插值是拉格朗日插值最簡單的情形)當(dāng)前8頁，總共42頁。分段三次埃爾米特插值條件數(shù)分段三次埃爾米特插值：線性插值在每一小段上(兩點(diǎn)之間)，用到2個(gè)條件q(xi)=yi，所以確定了一個(gè)線性插值函數(shù)；三次埃爾米特插值在每一小段上，用到4個(gè)條件q(xi)=yi，q'(xi)=y'i,所以確定一個(gè)3次多項(xiàng)式插值函數(shù)。分段插值主要是為了避免高次插值可能出現(xiàn)的大幅度振蕩現(xiàn)象，在實(shí)際應(yīng)用中通常采用分段低次插值來提高近似程度，比如可用分段線性插值或分段三次埃爾米特插值來逼近已知函數(shù)，但它們的總體光滑性較差，為了克服這一缺點(diǎn)，三次樣條插值成為比較理想的工具。當(dāng)前9頁，總共42頁。三次樣條(spline)插值的概念樣條的概念出自工程設(shè)計(jì)和機(jī)械加工(飛機(jī)、船舶外形曲線設(shè)計(jì))中的繪圖工具(曲線尺)，簡單說就是具有連續(xù)二階導(dǎo)數(shù)的三次插值多項(xiàng)式函數(shù)。

當(dāng)前10頁，總共42頁。三次樣條(spline)插值的條件數(shù)首先從段數(shù)n=2分析：我們知道在每一小段的三次多項(xiàng)式有4個(gè)系數(shù)，所以如下圖，總共需要有4*2=8個(gè)方程來確定；由q(xi)=yi可以確定2*2=4個(gè)方程，又由內(nèi)部節(jié)點(diǎn)q1'(xi)=q2'(xi)和q1''(xi)=q2''(xi)可以確定2*(2-1)=2個(gè)方程，看來剩下的8-(4+2)=2個(gè)方程只有靠外部給定(邊界條件)了q1q2x0x1x2當(dāng)前11頁，總共42頁。一維曲線等距插值函數(shù)interpinterp'ssyntaxOne-dimensionalrtimeslonger

datainterpolationy=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('OriginalPoints');subplot(1,2,2);stem(yi(1:120));title('InterpolatedPoints');當(dāng)前12頁，總共42頁。一維曲線等距插值函數(shù)interp1interp1'ssyntaxOne-dimensionaldatainterpolationyi=interp1(x,y,xi,method)'nearest'Nearestneighborinterpolation'linear'Linearinterpolation(default)'spline'Cubicsplineinterpolation'cubic'PiecewisecubicHermiteinterpolation題例在一天24小時(shí)內(nèi),從零點(diǎn)開始每間隔2小時(shí)測得的環(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');holdon;xi=15.1;yi=interp1(x,y,xi,'spline'),xi=0:1/3600:24;yi=interp1(x,y,xi,'spline');plot(xi,yi,'b-');當(dāng)前13頁，總共42頁。二維曲面等距插值函數(shù)interp2interp2'ssyntaxTwo-dimensionaldatainterpolationZI=interp2(X,Y,Z,XI,YI,method)'nearest'Nearestneighborinterpolation'linear'Bilinearinterpolation(default)'spline'Cubicsplineinterpolation'cubic'Bicubucinterpolation當(dāng)前14頁，總共42頁。二維曲面等距插值函數(shù)interp2動(dòng)畫展示：三維空間中的曲面等距格點(diǎn)當(dāng)前15頁，總共42頁。二維曲面等距插值函數(shù)interp2題例粗糙山頂曲面的平滑處理(等距情形)loadmountain.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,xii,yii,'spline');

%三次樣條插值figure;surf(xii,yii,zii)%繪制平滑處理后的山頂曲面holdon;[xx,yy]=meshgrid(x,y);plot3(xx,yy,z+0.1,'ob');當(dāng)前16頁，總共42頁。二維曲面等距插值函數(shù)interp2題例粗糙山頂曲面的平滑處理(等距情形)當(dāng)前17頁，總共42頁。二維曲面散亂插值函數(shù)griddatagriddata'ssyntaxDatainterpolation

forscatteredpoints

ZI=griddata(x,y,z,XI,YI)[XI,YI,ZI]=griddata(x,y,z,xi,yi)[...]=griddata(...,method)'linear'Triangle-basedlinearinterpolation'cubic'Triangle-basedcubic(default)'nearest'Nearestneighbor'v4'MATLAB4griddatamethodMATLAB二維插值函數(shù)griddata，可以將平面或曲面上的散亂點(diǎn)插值為規(guī)則網(wǎng)格當(dāng)前18頁，總共42頁。二維曲面散亂插值函數(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,'o');holdonti=-2:.25:2;[XI,YI]=meshgrid(ti,ti);ZI=griddata(x,y,z,XI,YI);mesh(XI,YI,ZI);當(dāng)前19頁，總共42頁。二維曲面散亂插值函數(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)holdonxlin=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);axistight;當(dāng)前20頁，總共42頁。南半球氣旋變化的可視圖形當(dāng)前21頁，總共42頁。山區(qū)地貌的可視化圖形當(dāng)前22頁，總共42頁。水塔用水量估計(jì)通用程序通用程序tbp69.m可近似計(jì)算時(shí)間段內(nèi)的用水量格式為：tbp69(ts,tf)其中ts為起點(diǎn)時(shí)間，tf為終點(diǎn)時(shí)間當(dāng)前23頁，總共42頁。實(shí)驗(yàn)一：水塔用水量估計(jì)

That’sall~3Q!

當(dāng)前24頁，總共42頁。第四講插值與擬合之?dāng)M合(下)內(nèi)容：擬合是離散函數(shù)逼近的重要方法，利用它可通過函數(shù)在有限個(gè)點(diǎn)處的取值狀況，擬合出近似替代函數(shù)，進(jìn)而估算出函數(shù)在其他點(diǎn)處的近似值。目的：學(xué)習(xí)擬合的基本思想和方法，掌握Matlab

的多項(xiàng)式/一般擬合函數(shù)/曲線擬合工具箱要求：掌握Matlab擬合函數(shù)，處理擬合應(yīng)用問題了解基于最小二乘法則擬合的基本思想掌握擬合函數(shù)polyfitlsqcurvefitcurvefit掌握cftool曲線擬合工具箱(多目標(biāo)函數(shù)多法則)當(dāng)前25頁，總共42頁。關(guān)于數(shù)據(jù)擬合的兩個(gè)要素...在工程實(shí)踐和科學(xué)計(jì)算中，用某種經(jīng)驗(yàn)函數(shù)解析式y(tǒng)=f(x)來近似刻畫采集數(shù)據(jù)(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ù)和擬合法則是擬合的兩個(gè)關(guān)鍵要素！當(dāng)前26頁，總共42頁。引例1化合物濃度隨時(shí)間變化的規(guī)律：與插值面臨的問題相似，我們被要求去求解或預(yù)測表格中沒有的因變量取值，與插值的解決思路不同，我們試圖獲得比較完備的解決方案：設(shè)計(jì)并求出離散數(shù)據(jù)點(diǎn)的近似替代函數(shù)，有了近似函數(shù)解析式，就可以進(jìn)一步代值計(jì)算或作圖分析。

化合物濃度隨時(shí)間變化的規(guī)律為了揭示濃度y與時(shí)間t之間呈現(xiàn)的函數(shù)規(guī)律，我們首先作出散點(diǎn)圖，幫助分析和設(shè)計(jì)經(jīng)驗(yàn)函數(shù)當(dāng)前27頁，總共42頁?；衔餄舛入S時(shí)間變化的規(guī)律如圖,化合物濃度y隨時(shí)間t大致呈拋物線狀（二次函數(shù)）變化，這種分析和判斷來自已有經(jīng)驗(yàn)...t=1:16;c=[46.488.49.289.59.79.861010.210.3210.4210.510.5510.5810.6];plot(t,c,'-ro')當(dāng)前28頁，總共42頁?；衔餄舛入S時(shí)間變化的規(guī)律經(jīng)驗(yàn)函數(shù)形式：已經(jīng)擬定為多項(xiàng)式函數(shù)：y=at2+bt+c剩下的工作是確定擬合原則：可選的法則很多，其中最常用的是最小二乘法則(methodofLeastSquares)，即各點(diǎn)殘差平方和最小高斯和勒讓德關(guān)于最小二乘法的發(fā)明權(quán)當(dāng)前29頁，總共42頁?；衔餄舛入S時(shí)間變化的規(guī)律對(duì)經(jīng)驗(yàn)函數(shù)形式確定的補(bǔ)充說明：擬合函數(shù)解析式選用什么形式(用多項(xiàng)式，還是用冪函數(shù)?)，主要取決于采樣點(diǎn)的分布無疑，那么如何求出這些含有待定參數(shù)的解析式呢？把各點(diǎn)偏差的平方和最小作為一個(gè)目標(biāo)函數(shù)，實(shí)際上考慮為極值問題，極值點(diǎn)導(dǎo)數(shù)為零。具體計(jì)算時(shí)，我們在把經(jīng)驗(yàn)函數(shù)用一系列擬合基函數(shù)線性表出同時(shí)，在j個(gè)采樣點(diǎn)對(duì)待定參數(shù)Cj求偏導(dǎo)(=0),獲取j個(gè)方程，進(jìn)而解出Cj，具體參見數(shù)值計(jì)算SZJSp90~91在本例中，已經(jīng)擬定擬合的目標(biāo)函數(shù)為多項(xiàng)式函數(shù)：y=at2+bt+c，所以只要解出三個(gè)待定參數(shù)a,b,c，問題即獲解決…當(dāng)前30頁，總共42頁?；谧钚《说亩囗?xiàng)式擬合函數(shù)基于最小二乘的多項(xiàng)式擬合函數(shù)polyfit：Polynomialcurvefitting.Syntax：p=polyfit(x,y,n)其中n是擬合多項(xiàng)式的階數(shù)，不能超過（散點(diǎn)數(shù)據(jù)對(duì)數(shù)-1）下面回到化合物濃度隨時(shí)間變化的引例：t=1:16;c=[46.488.49.289.59.79.861010.210.3210.4210.510.5510.5810.6];plot(t,c,'ko');holdon;%作散點(diǎn)圖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-');disp(['化合物在刻度11.2的濃度近似值為',num2str(f(11.2))])disp(['化合物在刻度17.8的濃度預(yù)測值為',num2str(f(17.8))])stem([11.217.8],[f(11.2)f(17.8)],'r');xlabel('時(shí)間t');ylabel('化合物濃度c');title('化合物濃度隨時(shí)間變化的規(guī)律')當(dāng)前31頁，總共42頁。引例2確定醫(yī)用薄膜滲透率的數(shù)學(xué)模型：

某種醫(yī)用薄膜允許一種物質(zhì)分子從高濃度溶液VB穿過薄膜向低濃度溶液VA中擴(kuò)散。通過單位面積膜S分子擴(kuò)散的速度與膜兩側(cè)溶液的濃度差成正比，比例系數(shù)K表示薄膜被該物質(zhì)分子穿透的能力，稱為滲透率，定時(shí)測量薄膜VB側(cè)的溶液濃度值CB，以此確定K的值VA=VB=1000cm3,

S=10cm2,

容器的B部分溶液濃度CB的測試結(jié)果如下表：(CB單位為mg/cm3

)確定醫(yī)用薄膜滲透率的數(shù)學(xué)模型當(dāng)前32頁，總共42頁。確定醫(yī)用薄膜滲透率的數(shù)學(xué)模型由質(zhì)量守恒考察[t,t+△t]時(shí)間段B向A中滲透物質(zhì)：VA*CA(t+△t)-VA*CA(t)=SK[CB(t)-CA(t)]△t

推出dCA(t)/dt=SK/VA*[CB(t)-CA(t)]

兩邊除以△t,令△t→0又由質(zhì)量守恒考察整個(gè)容器中物質(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當(dāng)前33頁，總共42頁。確定醫(yī)用薄膜滲透率的數(shù)學(xué)模型在上式中，代入已知值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.mfunctionCB=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é)果：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)當(dāng)前34頁，總共42頁?；谧钚《说囊话銛M合函數(shù)基于最小二乘的一般擬合函數(shù)lsqcurvefit：Solvenonlinearcurve-fitting(data-fitting)problemsintheleast-squaressense.Syntax：x=lsqcurvefit(fun,x0,xdata,ydata)[x,resnorm]=lsqcurvefit(...)范例:functionF=myfun(x,xdata)F=x(1)*xdata.^2+x(2)*sin(xdata)+x(3)*xdata.^3;下面是主調(diào)函數(shù)fitmyfun.m（片段）xdata=[3.67.79.34.18.62.81.37.910.05.4];ydata=[16.5150.6263.124.7208.59.92.7163.9325.054.3];x0=[10,10,10];%Startingguess[x,resnorm]=lsqcurvefit(@myfun,x0,xdata,ydata)當(dāng)前35頁，總共42頁?；谧钚《说囊话銛M合函數(shù)擬合:定參：x(1)=0.2269；x(2)=0.3385；x(3)=0.3021當(dāng)前36頁，總共42頁?；谧钚《说那鏀M合函數(shù)基于最小二乘的多元擬合函數(shù)lsqcurvefit：Solvenonlinearcurve-fitting(data-fitting)problemsintheleast-squaressense.Syntax：x=lsqcurvefit(fun,x0,xdata,ydata)[x,resnorm]=lsqcurvefit(...)范例:functionQ=mas(x,xdata)Q=x(1)*(xdata(1,:).^x(2)).*(xdata(2,:).^x(3));下面是主