數(shù)學(xué)實(shí)驗(yàn)8-曲線擬合與插值

MathematicsLaboratory阮小娥博士ExperimentsinMathematics數(shù)學(xué)實(shí)驗(yàn)在實(shí)際問(wèn)題中，我們常常會(huì)遇到下列問(wèn)題（1）變量間存在函數(shù)關(guān)系，但只能給出一離散點(diǎn)列上的值.例如，從實(shí)驗(yàn)中得到一個(gè)數(shù)據(jù)表,或是一組觀測(cè)數(shù)據(jù).（2）變量間的函數(shù)關(guān)系可以表示,但計(jì)算復(fù)雜,只能計(jì)算特殊點(diǎn)的函數(shù)值.例如，求指數(shù)函數(shù)、對(duì)數(shù)函數(shù)、三角函數(shù)、反三角函數(shù)值等.為了研究自變量與因變量間的變化關(guān)系，我們需要建立變量間的函數(shù)關(guān)系，從而可以計(jì)算原始數(shù)據(jù)以外需要處的值.解決這類問(wèn)題的方法：數(shù)據(jù)擬合、數(shù)據(jù)插值實(shí)驗(yàn)13人口數(shù)量預(yù)測(cè)模型實(shí)驗(yàn)2、掌握在最小二乘意義下數(shù)據(jù)擬合的理論和方法.1、學(xué)會(huì)用MATLAB軟件進(jìn)行數(shù)據(jù)擬合3、通過(guò)對(duì)實(shí)際問(wèn)題的分析和研究，初步掌握建立數(shù)據(jù)擬合數(shù)學(xué)模型的方法實(shí)驗(yàn)?zāi)康膿?jù)人口統(tǒng)計(jì)年鑒，知我國(guó)從1949年至1994年人口數(shù)據(jù)資料如下：(人口數(shù)單位為：百萬(wàn))（1）在直角坐標(biāo)系上作出人口數(shù)的圖象。（2）建立人口數(shù)與年份的函數(shù)關(guān)系，并估算1999年的人口數(shù)。實(shí)驗(yàn)問(wèn)題年份

1949

1954

1959

1964

1969

人口數(shù)

541.67

602.66

672.09704.99

806.71

年份1974

1979

1984

1989

1994人口數(shù)

908.59

975.421034.75

1106.76

1176.74

如何確定a,b?線性模型1曲線擬合問(wèn)題的提法:

已知一組（二維）數(shù)據(jù)，即平面上的n個(gè)點(diǎn)),(iiyx，

ixni,,,2,1L=互不相同，尋求一個(gè)函數(shù)（曲線）)(xfy=，使)(xf在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好，如圖:

xy0++++++++一、曲線擬合確定f(x)使得

達(dá)到最小

最小二乘準(zhǔn)則

２.用什么樣的曲線擬合已知數(shù)據(jù)?常用的曲線函數(shù)系類型：(１)畫(huà)圖觀察(２)理論分析指數(shù)曲線：

雙曲線（一支):

多項(xiàng)式：

直線：

３擬合函數(shù)組中系數(shù)的確定二、人口預(yù)測(cè)的線性模型

對(duì)于開(kāi)始提出的實(shí)驗(yàn)問(wèn)題,代如數(shù)據(jù)，計(jì)算得從而得到人口數(shù)與年份的函數(shù)關(guān)系為把x=1999代如，估算出1999年的人口數(shù)為y=1252.1（百萬(wàn)）＝12.52億1999年實(shí)際人口數(shù)量為１２.６億。線性預(yù)測(cè)模型英國(guó)人口學(xué)家Malthus根據(jù)百余年的人口統(tǒng)計(jì)資料，于1798年提出了著名的人口自然增長(zhǎng)的指數(shù)增長(zhǎng)模型。三、人口預(yù)測(cè)的Malthus模型基本假設(shè)

:人口(相對(duì))增長(zhǎng)率r

是常數(shù)x(t)~時(shí)刻t的人口,t=0時(shí)人口數(shù)為x0指數(shù)增長(zhǎng)模型實(shí)際中，常用1.由前100年的數(shù)據(jù)求出美國(guó)的人口增長(zhǎng)Malthus模型。2.預(yù)測(cè)后100年（每隔10年）的人口狀況。3.根據(jù)預(yù)測(cè)的人口狀況和實(shí)際的人口數(shù)量,討論人口模型的改進(jìn)情況。美國(guó)1790年－1980年每隔10年的人口記錄226.5204.0179.3150.7131.7123.2106.592.076.062.9人口(百萬(wàn))1980197019601950194019301920191019001890年份50.238.631.423.217.112.99.67.25.33.9人口(百萬(wàn))1880187018601850184018301820181018001790年份例１解：取得最小值.其中,表示人口數(shù)量。表示年份,解方程組:即得參數(shù)的值.使得問(wèn)題轉(zhuǎn)化為求參數(shù)％prog41.m%%Thisprogramistopredictthenumberofpopulation%formatlongt1=[1790;1800;1810;1820;1830;1840;1850;1860;1870;1880];t2=[1890;1900;1910;1920;1930;1940;1950;1960;1970;1980];x1=[3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2];x2=[62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5];lnx1=log(x1);lnx2=log(x2);a12=sum(t1);a11=10;a21=a12;a22=sum(t1.^2);d1=sum(lnx1);d2=sum(lnx1.*t1);

A=[a11,a12;a21,a22];D=[d1;d2];

ab=inv(A)*D;

disp('a=');disp(ab(1));

disp('b=');disp(ab(2));

fori=1:10

xx1(i)=exp(ab(1)+ab(2)*t1(i));

end

fori=1:10

xx2(i)=exp(ab(1)+ab(2)*t2(i));

end

plot(t1,x1,'r*--',t1,xx1,'b+-',t2,x2,'g*--',t2,xx2,'m+-');

a=-49.79535457790735b=0.02859807120038仿真結(jié)果表明：人口增加的指數(shù)模型在短期內(nèi)基本上能比較準(zhǔn)確地反映人口自然增長(zhǎng)的規(guī)律，但長(zhǎng)期預(yù)測(cè)誤差很大，需要修正預(yù)測(cè)模型。擬合曲線原始數(shù)據(jù)曲線四、人口預(yù)測(cè)的Logistic模型人口增長(zhǎng)到一定數(shù)量后，增長(zhǎng)率下降的原因：資源、環(huán)境等因素對(duì)人口增長(zhǎng)的阻滯作用且阻滯作用隨人口數(shù)量增加而變大假設(shè)r~固有增長(zhǎng)率(x很小時(shí))k~人口容量（資源、環(huán)境能容納的最大數(shù)量）r是x的減函數(shù)例１的Logistic模型留給同學(xué)們練習(xí)五、多項(xiàng)式擬合的Matlab指令a=polyfit（xdata，ydata，n）其中n表示多項(xiàng)式的最高階數(shù)

xdata，ydata為要擬合的數(shù)據(jù)，它是用向量的方式輸入。輸出參數(shù)a為擬合多項(xiàng)式

y=a1xn+…+anx+an+1的系數(shù)a=[a1,…,an,an+1]。多項(xiàng)式在x處的值y可用下面程序計(jì)算。

y=polyval(a,x)

用多項(xiàng)式擬合人口模型%Thisprogramistopredictthemodelofpopulationby4-degreepolynomial%%prog42.m%formatlongt1=[1790;1800;1810;1820;1830;1840;1850;1860;1870;1880];t2=[1890;1900;1910;1920;1930;1940;1950;1960;1970;1980];t=[t1;t2];P1=[3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2];P2=[62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5];P=[P1;P2];n=4;%Thedegreeofthefittingpolynomial%[a,s]=polyfit(t1,P1,n);y=polyval(a,t);%aisthecoefficientsvectorfromn-degreeto0-degree%plot(t,P,'r*--',t,y,'b+-');23a=1.0e+006*-0.000000000000140.00000000107892-0.000003048785950.00381927346813-1.79012132225427仿真結(jié)果表明,人口增加的模型用多項(xiàng)式擬合能比較準(zhǔn)確地反映人口自然增長(zhǎng)的規(guī)律，對(duì)長(zhǎng)期預(yù)測(cè)具有指導(dǎo)意義。例2:海底光纜線長(zhǎng)度預(yù)測(cè)模型某一通信公司在一次施工中,需要在水面寬為20m的河溝底沿直線走向鋪設(shè)一條溝底光纜.在鋪設(shè)光纜之前需要對(duì)溝底的地形做初B2468101214161820986420ADC探測(cè)到一組等分點(diǎn)位置的深度數(shù)據(jù)如下表所示.25步探測(cè),從而估計(jì)所需光纜的長(zhǎng)度,為工程預(yù)算提供依據(jù).基本情況如圖所示.10.9310.809.818.867.957.959.1510.2211.2912.6113.32201918171615141312111013.2812.2611.1810.139.058.027.967.968.969.01深度(m)9876543210分點(diǎn)21個(gè)等分點(diǎn)處的深度(1)預(yù)測(cè)通過(guò)這條河溝所需光纜長(zhǎng)度的近似值.(2)作出鋪設(shè)溝底光纜的曲線圖.解：用12次多項(xiàng)式函數(shù)擬合光纜走勢(shì)的曲線圖如下仿真結(jié)果表明,擬合曲線能較準(zhǔn)確地反映光纜的走勢(shì)圖.Thelengthofthelabelis

L=26.3809(m)假設(shè)所鋪設(shè)的光纜足夠柔軟,在鋪設(shè)過(guò)程中光纜觸地走勢(shì)光滑,緊貼地面,并且忽略水流對(duì)光纜的沖擊.%prog45.mThisprogramistofitthedatabypolynomial%formatlongt=linspace(0,20,21);x=linspace(0,20,100);P=[9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93];[a,s]=polyfit(t,P,12);yy=polyval(a,x);plot(x,yy,'r*--',t,P,'b+-');L=0;fori=2:100L=L+sqrt((x(i)-x(i-1))^2+(yy(i)-yy(i-1))^2);enddisp('ThelengthofthelabelisL=');disp(L);formatlongt=linspace(0,20,21);x=linspace(0,20,100);P=[9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93];n=input(‘n=’)%通過(guò)鍵盤(pán)輸入擬合次數(shù)[a,s]=polyfit(t,P,n);yy=polyval(a,x);p1=polyval(a,t);d=norm(P-p1)%計(jì)算擬合誤差plot(x,yy,'r*--',t,P,'b+-');L=0;fori=2:100L=L+sqrt((x(i)-x(i-1))^2+(yy(i)-yy(i-1))^2);enddisp('ThelengthofthelabelisL=');disp(L);六、最小二乘曲線擬合有一組數(shù)據(jù)xi，yi（i=1,2,…,n）滿足某一函數(shù)原型，其中a為待定系數(shù)向量，求出一組待定系數(shù)的值使得目標(biāo)函數(shù)最?。鹤钚《饲€擬合函數(shù)lsqcurvefit的調(diào)用格式：[a,Jm]=lsqcurvefit(Fun,a0,x,y)Fun為原型函數(shù)的matlab表示，可以是M-函數(shù)或inline()函數(shù)a0為最優(yōu)化初值x和y為原始輸入輸出數(shù)據(jù)向量a為返回的待定系數(shù)向量Jm為在待定系數(shù)下目標(biāo)函數(shù)的值例3

已知數(shù)據(jù)可能滿足求滿足數(shù)據(jù)的最小二乘解a、b、c和d的值.x0.10.20.30.40.5y2.32012.64702.97073.28853.6008x0.60.70.80.91.0y3.90904.21474.51914.82325.1275輸入已知數(shù)據(jù)點(diǎn)：x=0.1:

0.1

:

1;y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,4.2147,4.5191,4.8232,5.1275];編寫(xiě)函數(shù)functiony=f3(a,x)y=a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4);待定系數(shù)求解a=lsqcurvefit('f3',[1;2;2;3],x,y);a'ans=2.45752.45571.44372.0720繪制曲線：>>y1=f3(a,x);plot(x,y,x,y1,’o’)插值問(wèn)題：實(shí)驗(yàn)14插值問(wèn)題插值條件插值函數(shù)插值節(jié)點(diǎn)如果是多項(xiàng)式,則稱為插值多項(xiàng)式求插值函數(shù)的方法稱為插法…………[a,b]稱為插值區(qū)間如何構(gòu)造P(x)?1、多項(xiàng)式插值法當(dāng)n=0時(shí),只有一個(gè)插值節(jié)點(diǎn)的情形當(dāng)n=1時(shí),有兩個(gè)插值節(jié)點(diǎn)的情形當(dāng)n=2時(shí),有三個(gè)插值節(jié)點(diǎn)的情形是否任意給定n+1個(gè)不同的插值節(jié)點(diǎn)都可以構(gòu)造出滿足插值條件的插值多項(xiàng)式?由克萊姆法則知方程組有唯一解，即滿足(1.1)的插值多項(xiàng)式存在且唯一。插值多項(xiàng)式的存在唯一定理：在次數(shù)不超過(guò)n的多項(xiàng)式集合中，滿足插值條件的插值多項(xiàng)式是存在并且唯一的。拉格朗日(Lagrange)插值多項(xiàng)式表示插值多項(xiàng)式的最緊湊的方式是拉格朗日形式：

Lagrange插值多項(xiàng)式的優(yōu)點(diǎn)在于不要求數(shù)據(jù)點(diǎn)是等間隔的，其缺點(diǎn)是數(shù)據(jù)點(diǎn)數(shù)不宜過(guò)大，通常不超過(guò)7個(gè)，否則計(jì)算工作量大且誤差大，計(jì)算不穩(wěn)定。分段線性插值分段線性插值示意圖…………分段二次插值示意圖…………分段二次插值三次樣條插值函數(shù)定義對(duì)于區(qū)間[a,b]上給定的一個(gè)分劃如果函數(shù)S(x)在子區(qū)間上都是不超過(guò)3次的多項(xiàng)式，并且2階導(dǎo)數(shù)在內(nèi)節(jié)點(diǎn)處連續(xù)，則稱為區(qū)間[a,b]上以為節(jié)點(diǎn)的三次樣條函數(shù)。對(duì)于函數(shù)，若還滿足插值條件：則稱為在區(qū)間上的三次樣條插值函數(shù)。是在飛機(jī)或輪船等的設(shè)計(jì)制造過(guò)程中為描繪出光滑的外形曲線(放樣)所用的工具.2、插值函數(shù)的MATLAB實(shí)現(xiàn)（1）interp1函數(shù)interp1函數(shù)的調(diào)用格式：y=interp1(x0,y0,x,method)其中:1）x0、y0為樣本點(diǎn)，y為插值點(diǎn)自變量值x對(duì)應(yīng)的函數(shù)值。（2）method共有6種參數(shù)可供選擇，當(dāng)省略method時(shí)，即默認(rèn)為linear線性插值。線性插值：method=‘linear’三次樣條插值:method=‘spline’;立方插值：method=‘pchip’or‘cubic’例4已知某轉(zhuǎn)子流量計(jì)在100～1000mL/min流量范圍內(nèi),刻度值與校正值的關(guān)系如表所示.試用線性插值法計(jì)算流量計(jì)的刻度值為785時(shí),實(shí)際流量為多少？刻度值校正值刻度值校正值100105.3600605.8200207.2700707.4300308.1800806.7400406.9900908.0500507.51000107.9解：Matlab計(jì)算程序如下：X=[100,200,300,400,500,600,700,800,900,1000];Y=[105.3,207.2,308.1,406.9,507.5,605.8,707.4,806.7,908.0,107.9];Xk=780;Yk=interp1(X,Y,Xk)執(zhí)行結(jié)果：Yk=786.8400這里：X和Y分別表示樣本點(diǎn)的刻度值和校正值；

Xk和Yk分別表示插值點(diǎn)的刻度值和校正值。功能

三次樣條數(shù)據(jù)插值格式

y=spline(x0,y0,x)

與y=interp1(x0,y0,x,’spline’)等價(jià)，其中參數(shù)x0、y0、x，y的意義及要求與線性插值interp1中的完全一致。（2）spline函數(shù)例5對(duì)離散分布在y=exp(x)sin(x)函數(shù)曲線上的數(shù)據(jù)點(diǎn)進(jìn)行樣條插值計(jì)算：

x=[024581212.817.219.920];y=exp(x).*sin(x);xx=0:0.25:20;yy=spline(x,y,xx);plot(x,y,'bo',xx,yy,'r*')例6

已知某型號(hào)飛機(jī)的機(jī)翼斷面下緣輪廓線上的部分?jǐn)?shù)據(jù)如表所示:假設(shè)需要得到x

坐標(biāo)每改變0.1時(shí)的y

坐標(biāo),

分別用兩種插值方法對(duì)機(jī)翼斷面下緣輪廓線上的部分?jǐn)?shù)據(jù)加細(xì),并作出插值函數(shù)的圖形.xy0031.251.772.092.1112.0121.8131.2141.0151.6程序如下：x=[0,3,5,7,9,11,12,13,14,15]y=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6]xi=[0:0.1:15]yi=interp1(x,y,xi,'spline')zi=spline(x,y,xi)plot(xi,yi,’bO’,xi,zi,’r*’)（3）lagrange插值函數(shù)function

y=lagrange(x0,y0,x)ii=1:length(x0);y=zeros(size(x));fori=iiij=find(ii~=i);y1=1;forj=1:length(ij),y1=y1.*(x-x0(ij(j)));endy=y+y1*y0(i)/prod(x0(i)-x0(ij));end例7對(duì)進(jìn)行Lagrange插值x0=-1+2*[0:10]/10;y0=1./(1+25*x0.^2);x=-1:0.01:1;y=lagrange(x0,y0,x);ya=1./(1+25*x.^2);

plot(x,ya,'*',x,y,'O')3、二維網(wǎng)格數(shù)據(jù)的插值問(wèn)題（1）二維插值函數(shù)interp2的調(diào)用格式：zi=interp2(x0,y0,z0,xi,yi)zi=interp2(x0,y0,z0,xi,yi,

method)*臨近點(diǎn)插值：method=‘nearest’*線性插值：method=‘linear’（缺省算法）*三次樣條插值：method=‘spline’*立方插值：method=‘pchip’or‘cubic’例8

由二元函數(shù)獲得一些較稀疏的網(wǎng)格數(shù)據(jù)，對(duì)整個(gè)函數(shù)曲面進(jìn)行各種插值，并比較插值結(jié)果%繪制已知數(shù)據(jù)的網(wǎng)格圖[x,y]=meshgrid(-3:0.6:3,-2:0.4:2);z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);surf(x,y,z);axis([-3,3,-2,2,-0.7,1.5])%默認(rèn)線性插值算法進(jìn)行插值[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);z1=interp2(x,y,z,x1,y1);surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])%立方和樣條插值z(mì)2=interp2(x,y,z,x1,y1,'cubic');z3=interp2(x,y,z,x1,y1,'spline');surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])figure;surf(x1,y1,z3),axis([-3,3,-2,2,-0.7,1.5])[x,y]=meshgrid(-3:0.6:3,-2:0.4:2);z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);subplot(2,2,1);surf(x,y,z);axis([-3,3,-2,2,-0.7,1.5])[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);z1=interp2(x,y,z,x1,y1);subplot(2,2,2);surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])z2=interp2(x,y,z,x1,y1,'cubic');z3=interp2(x,y,z,x1,y1,'spline');subplot(2,2,3);surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])subplot(2,2,4);surf(x1,y1,z3),axis([-3,3,-2,2,-0.7,1.5])畫(huà)在同一幅圖上作比較

前兩個(gè)差值算法誤差的比較z=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);surf(x1,y1,abs(z-z1)),axis([-3,3,-2,2,0,0.08])figure;surf(x1,y1,abs(z-z2)),axis([-3,3,-2,2,0,0.025])三個(gè)差值算法誤差的比較z=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);subplot(3,1,1);surf(x1,y1,abs(z-z1)),axis([-3,3,-2,2,0,0.18])subplot(3,1,2);surf(x1,y1,abs(z-z2)),axis([-3,3,-2,2,0,0.18])subplot(3,1,3);surf(x1,y1,abs(z-z3)),axis([-3,3,-2,2,0,0.18])（2）二維一般分布數(shù)據(jù)的插值問(wèn)題griddata函數(shù)的調(diào)用格式：z=griddata(x0,y0,z0,x,y,method)

method=‘v4’：插值算法，公認(rèn)效果較好臨近點(diǎn)插值：method=‘nearest’線性插值：method=‘linear’（缺省算法）三次樣條插值：method=‘spline’立方插值：method=‘cubic’例9

在x為［-3,3]，y為［-2，2］矩形區(qū)域隨機(jī)選擇一組數(shù)據(jù)點(diǎn)，用’v4’與’cubic’插值法進(jìn)行處理，并對(duì)誤差進(jìn)行比較。%產(chǎn)生數(shù)據(jù)點(diǎn)x=-3+6*rand(200,1);y=-2+4*rand(200,1);z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);plot(x,y,'x')%樣本點(diǎn)的二維分布

figure,plot3(x,y,z,'x'),axis([-3,3,-2,2,-0.7,1.5]),grid%cubic和v4算法[x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2);z1=griddata(x,y,z,x1,y1,'cubic');subplot(2,2,1);surf(x1,y1,z1);axis([-3,3,-2,2,-0.7,1.5])z2=griddata(x,y,z,x1,y1,'v4');subplot(2,2,2);surf(x1,y1,z2);axis([-3,3,-2,2,-0.7,1.5])%誤差分析z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);subplot(2,2,3);surf(x1,y1,abs(z0-z1)),axis([-3,3,-2,2,0,0.15])subplot(2,2,4);surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.15])4、高維插值問(wèn)題三維插值interp3函數(shù)的調(diào)用格式：三維網(wǎng)格[x,y,z]=meshgrid(x1,y1,z1)griddata3()三維非網(wǎng)格數(shù)據(jù)插值n維插值interpn函數(shù)n維網(wǎng)格[x1,x2,…,xn]=ndgrid[v1,v2,…,vn]griddatan()n維非網(wǎng)格數(shù)據(jù)插值interp3()、interpn()調(diào)用格式同interp2()一致griddata3()、griddatan()調(diào)用格式同griddata()一致例10

通過(guò)函數(shù)生成一些網(wǎng)格型樣本點(diǎn)，據(jù)此進(jìn)行插值并給出插值誤差。[x,y,z]=meshgrid(-1:0