matlab人口預(yù)測(cè)與數(shù)據(jù)曲線擬合

MathematicsLaboratoryExperimentsinMathematics趙小艷數(shù)學(xué)實(shí)驗(yàn)辦公地址：理科樓2142021/5/91實(shí)驗(yàn)13人口預(yù)測(cè)與數(shù)據(jù)擬合2、了解利用最小二乘法進(jìn)行數(shù)據(jù)擬合的基本思想，掌握用數(shù)據(jù)擬合法尋找最佳擬合曲線的方法。3、了解多元函數(shù)的極值在數(shù)據(jù)擬合法中的應(yīng)用。實(shí)驗(yàn)?zāi)康?、學(xué)會(huì)用MATLAB軟件進(jìn)行數(shù)據(jù)擬合。4、通過對(duì)實(shí)際問題進(jìn)行分析研究，初步掌握建立數(shù)據(jù)擬合數(shù)學(xué)模型的方法。2021/5/92據(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)問題年份

19491954195919641969

人口數(shù)

541.67602.66672.09704.99806.71

年份19741979198419891994人口數(shù)

908.59975.421034.751106.761176.74

2021/5/93如何確定a,b?線性模型2021/5/941曲線擬合問題的提法:

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

ixni,,,2,1L=互不相同，尋求一個(gè)函數(shù)（曲線）)(xfy=，使)(xf在觀測(cè)點(diǎn)x處所取得值f(x)分別與觀察值y在某種

xy0++++++++一、曲線擬合準(zhǔn)則下最為接近，即曲線擬合得最好，如圖2021/5/95從幾何上講，并不要求曲線嚴(yán)格通過已知點(diǎn)，但要求曲線在各數(shù)據(jù)點(diǎn)和已知數(shù)據(jù)點(diǎn)之間的總體誤差最小，通常稱為數(shù)據(jù)擬合。達(dá)到最小。

最小二乘準(zhǔn)則

而我們經(jīng)常是確定f(x)使得偏差平方和2021/5/96數(shù)據(jù)插值已知一組（二維）數(shù)據(jù)，即平面上的n個(gè)點(diǎn)),(iiyx，

ixni,,,2,1L=互不相同，尋求一個(gè)函數(shù)（曲線）)(xfy=

數(shù)據(jù)插值2021/5/97２.用什么樣的曲線擬合已知數(shù)據(jù)?常用的曲線函數(shù)系ri(x)類型：１）畫圖觀察２）理論分析指數(shù)曲線：

雙曲線（一支):

多項(xiàng)式：

直線：

2021/5/98例如指數(shù)函數(shù)擬合三角函數(shù)擬合多項(xiàng)式擬合2021/5/99３擬合函數(shù)組中系數(shù)的確定2021/5/9104用matlab軟件進(jìn)行數(shù)據(jù)擬合（1）lsqcurvefit命令最小二乘擬合a=lsqcurvefit(fun,x0,xdata,ydata)[a,resnorm]=lsqcurvefit(fun,x0,xdata,ydata)

是根據(jù)給定的數(shù)據(jù)xdata,ydata，按照函數(shù)文件fun給定的函數(shù)，以x0為初值做最小二乘擬合，返回函數(shù)中的系數(shù)向量a和殘差平方和resnorm。2021/5/911例首先編寫函數(shù)文件functiony=f(a,x)f=a(1)*exp(x)+a(2)*x.^2+a(3)*x.^3保存為f.m，其次調(diào)用該函數(shù)x=0:0.1:1;y=[3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17];a0=[000];[x,resnorm]=lsqcurvefit(@f,a0,x,y)2021/5/912也可以用inline命令定義函數(shù)x=0:0.1:1;y=[3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17];f=inline('a(1)*exp(x)+a(2)*x.^2+a(3)*x.^3','a','x');a0=[000];[a,resnorm]=lsqcurvefit(f,a0,x,y)plot(x,y,'*')holdong=a(1)*exp(x)+a(2)*x.^2+a(3)*x.^3;plot(x,g,'r-')2021/5/913a=polyfit(xdata,ydata,n)其中n表示多項(xiàng)式的最高階數(shù)

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

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

y=polyval(a,x)

由于高次多項(xiàng)式曲線變化不穩(wěn)定，所以多項(xiàng)式次數(shù)的選取不宜過高。(2)polyfit命令多項(xiàng)式曲線擬合2021/5/914例如clear;clc;x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.66,9.56,9.48,9.3,11.2];plot(x,y,'k.','markersize',25);axis([01.3-216]);p3=polyfit(x,y,3)p6=polyfit(x,y,6)t=0:0.01:1.2;s=polyval(p3,t);s1=polyval(p6,t);holdonplot(t,s,'r-','linewidth',2);plot(t,s1,'b--','linewidth',2);grid2021/5/9152021/5/916編寫程序調(diào)用matlab命令x=1949:5:1994;y=[541.67,602.66,672.09,704.99,806.71,908.59,975.42,1034.75,1106.76,1176.74];plot(x,y,'r*','linewidth',2)gridf=inline('a(1)+a(2)*x','a','x');a0=[05];[a,resnorm]=lsqcurvefit(f,a0,x,y)holdong=a(1)+a(2)*x;plot(x,g,'b-','linewidth',2)二、人口預(yù)測(cè)線性模型2021/5/917或者調(diào)用M函數(shù)functionf=nihe(a,x)f=a(1)+a(2)*x;保存成nihe.m,在新窗口編寫程序x=1949:5:1994;y=[541.67,602.66,672.09,704.99,806.71,908.59,975.42,1034.75,1106.76,1176.74];a0=[1010];[a,resnorm]=lsqcurvefit(@nihe,a0,x,y)注意：該命令與初值有關(guān)系。2021/5/918也可以直接編寫程序如下：clc;clf;x=1949:5:1994;y=[541.67,602.66,672.09,704.99,806.71,908.59,975.42,1034.75,1106.76,1176.74];plot(x,y,'r*','linewidth',2)grida11=10;a12=sum(x);a21=a12;a22=sum(x.^2);d1=sum(y);d2=sum(x.*y);A=[a11,a12;a21,a22];D=[d1;d2];ab=inv(A)*Dplot(x,g,'b-','linewidth',2)t=1949:5:2010;g=ab(1)+ab(2)*t;holdonplot(t,g,'b-','linewidth',2)y2000=ab(1)+ab(2)*2000y2005=ab(1)+ab(2)*2005y2010=ab(1)+ab(2)*2010axis([194520125001450])plot(2000,1295.3,'g*','linewidth',2)plot(2005,1306.28,'g*','linewidth',2)plot(2010,1370.5,'g*','linewidth',2)2021/5/919計(jì)算得從而得到人口數(shù)與年份的函數(shù)關(guān)系為線性預(yù)測(cè)模型年份200020052010預(yù)測(cè)（百萬(wàn)）1266.61339.11411.7真實(shí)值（百萬(wàn)）1295.331306.281370.5并預(yù)測(cè)2000,2005,2010年的人口2021/5/920仿真結(jié)果表明：線性模型在短期內(nèi)基本上能比較準(zhǔn)確地反映人口自然增長(zhǎng)的規(guī)律，但長(zhǎng)期預(yù)測(cè)誤差較大。2021/5/921英國(guó)統(tǒng)計(jì)學(xué)家Malthus于1798年提出了一種關(guān)于生物種群繁殖的指數(shù)增長(zhǎng)模型：假設(shè)種群數(shù)量的增長(zhǎng)率與該時(shí)刻種群的個(gè)體數(shù)量成正比。三、人口預(yù)測(cè)的Malthus模型基本假設(shè)

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

是常數(shù)x(t)~時(shí)刻t的人口,t=0時(shí)人口數(shù)為x0指數(shù)增長(zhǎng)模型實(shí)際中，常用2021/5/922解：取得最小值.其中,表示人口數(shù)量。表示年份,解方程組:即得參數(shù)的值.使得問題轉(zhuǎn)化為求參數(shù)2021/5/923計(jì)算得從而得到人口數(shù)與年份的函數(shù)關(guān)系為指數(shù)預(yù)測(cè)模型年份200020052010預(yù)測(cè)（百萬(wàn)）1363.61488.81625.4真實(shí)值（百萬(wàn)）1295.331306.281370.5并預(yù)測(cè)2000,2005,2010年的人口2021/5/924仿真結(jié)果表明：人口增加的指數(shù)模型在較短期內(nèi)基本上能比較準(zhǔn)確地反映人口自然增長(zhǎng)的規(guī)律，但長(zhǎng)期預(yù)測(cè)誤差很大。2021/5/925四、人口預(yù)測(cè)的Logistic模型如果人口的增長(zhǎng)符合Malthus模型，則當(dāng)人口增長(zhǎng)到一定數(shù)量后，增長(zhǎng)率下降的原因：資源、環(huán)境等因素對(duì)人口增長(zhǎng)的阻滯作用1838年，荷蘭生物學(xué)家Verhulst對(duì)Malthus模型作了進(jìn)一步分析后指出：導(dǎo)致上述不符合實(shí)際情況的主要原因是未能考慮“密度制約”因素。即最終導(dǎo)致地球上人口爆炸，這與實(shí)際是不相符的。且阻滯作用隨人口數(shù)量增加而變大r是x的減函數(shù)2021/5/926假設(shè)r~固有增長(zhǎng)率(x很小時(shí))k~人口容量（資源、環(huán)境能容納的最大數(shù)量）2021/5/927結(jié)果顯示：Logistic模型在長(zhǎng)期預(yù)測(cè)時(shí)基本上能比較準(zhǔn)確地反映人口自然增長(zhǎng)的規(guī)律，在2005,2010年較為吻合。2021/5/928五、人口預(yù)測(cè)的多項(xiàng)式模型仿真結(jié)果表明,人口模型用低階多項(xiàng)式擬合能比較準(zhǔn)確地反映人口自然增長(zhǎng)的規(guī)律，而高階多項(xiàng)式擬合預(yù)測(cè)效果很差。青色二次多項(xiàng)式藍(lán)色—三次多項(xiàng)式紫紅色四次多項(xiàng)式2021/5/929例2:海底光纜線長(zhǎng)度預(yù)測(cè)模型某一通信公司在一次施工中,需要在水面寬為20m的河溝底沿直線走向鋪設(shè)一條溝底光纜.在鋪設(shè)光纜之前需要對(duì)溝底的地形做初B2468101214161820986420ADC探測(cè)到一組等分點(diǎn)位置的深度數(shù)據(jù)如下表所示.25步探測(cè),從而估計(jì)所需光纜的長(zhǎng)度,為工程預(yù)算提供依據(jù).基本情況如圖所示.2021/5/93010.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è)通過這條河溝所需光纜長(zhǎng)度的近似值.(2)作出鋪設(shè)溝底光纜的曲線圖.2021/5/931解：用12次多項(xiàng)式函數(shù)擬合光纜走勢(shì)的曲線圖如下仿真結(jié)果表明,擬合曲線能較準(zhǔn)確地反映光纜的走勢(shì)圖.ThelengthofthelabelisL=26.3809(m)假設(shè)所鋪設(shè)的光纜足夠柔軟,在鋪設(shè)過程中光纜觸地走勢(shì)光滑,緊貼地面,并且忽略水流對(duì)光纜的沖擊.2021/5/932%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);disp('yy=');disp(yy);plot(x,yy,