(五)-MATLAB應(yīng)用舉例資料課件

MATLAB應(yīng)用舉例例12005高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽題目C雨量預(yù)報(bào)方法的評(píng)價(jià)

雨量預(yù)報(bào)對(duì)農(nóng)業(yè)生產(chǎn)和城市工作和生活有重要作用，但準(zhǔn)確、及時(shí)地對(duì)雨量作出預(yù)報(bào)是一個(gè)十分困難的問(wèn)題，廣受世界各國(guó)關(guān)注。我國(guó)某地氣象臺(tái)和氣象研究所正在研究6小時(shí)雨量預(yù)報(bào)方法，即每天晚上20點(diǎn)預(yù)報(bào)從21點(diǎn)開(kāi)始的4個(gè)時(shí)段（21點(diǎn)至次日3點(diǎn)，次日3點(diǎn)至9點(diǎn)，9點(diǎn)至15點(diǎn)，15點(diǎn)至21點(diǎn)）在某些位置的雨量，這些位置位于東經(jīng)120度、北緯32度附近的53×47的等距網(wǎng)格點(diǎn)上。同時(shí)設(shè)立91個(gè)觀(guān)測(cè)站點(diǎn)實(shí)測(cè)這些時(shí)段的實(shí)際雨量，由于各種條件的限制，站點(diǎn)的設(shè)置是不均勻的。10/6/20230MATLAB應(yīng)用舉例例12005高教社杯全國(guó)大學(xué)生數(shù)學(xué)例1雨量預(yù)報(bào)方法的評(píng)價(jià)氣象部門(mén)希望建立一種科學(xué)評(píng)價(jià)預(yù)報(bào)方法好壞的數(shù)學(xué)模型與方法。氣象部門(mén)提供了41天的用兩種不同方法的預(yù)報(bào)數(shù)據(jù)和相應(yīng)的實(shí)測(cè)數(shù)據(jù)。預(yù)報(bào)數(shù)據(jù)在文件夾FORECAST中，實(shí)測(cè)數(shù)據(jù)在文件夾MEASURING中，其中的文件都可以用Windows系統(tǒng)的“寫(xiě)字板”程序打開(kāi)閱讀。FORECAST中的文件lon.dat和lat.dat分別包含網(wǎng)格點(diǎn)的經(jīng)緯度，其余文件名為<f日期i>_dis1和<f日期i>_dis2，例如f6181_dis1中包含2002年6月18日晚上20點(diǎn)采用第一種方法預(yù)報(bào)的第一時(shí)段數(shù)據(jù)（其2491個(gè)數(shù)據(jù)為該時(shí)段各網(wǎng)格點(diǎn)的雨量），而f6183_dis2中包含2002年6月18日晚上20點(diǎn)采用第二種方法預(yù)報(bào)的第三時(shí)段數(shù)據(jù)。10/6/20231例1雨量預(yù)報(bào)方法的評(píng)價(jià)氣象部門(mén)希望建立一種科學(xué)評(píng)價(jià)預(yù)報(bào)方法例1雨量預(yù)報(bào)方法的評(píng)價(jià)MEASURING中包含了41個(gè)名為<日期>.SIX的文件，如020618.SIX表示2002年6月18日晚上21點(diǎn)開(kāi)始的連續(xù)4個(gè)時(shí)段各站點(diǎn)的實(shí)測(cè)數(shù)據(jù)（雨量），這些文件的數(shù)據(jù)格式是：站號(hào)緯度經(jīng)度第1段第2段第3段第4段5813832.9833118.51670.00000.200010.10003.10005813933.3000118.85000.00000.00004.60007.40005814133.6667119.26670.00000.00001.10001.40005814333.8000119.80000.00000.00000.00001.80005814633.4833119.81670.00000.00001.50001.9000……雨量用毫米做單位，小于0.1毫米視為無(wú)雨。10/6/20232例1雨量預(yù)報(bào)方法的評(píng)價(jià)MEASURING中包例5.1雨量預(yù)報(bào)方法的評(píng)價(jià)（1）請(qǐng)建立數(shù)學(xué)模型來(lái)評(píng)價(jià)兩種6小時(shí)雨量預(yù)報(bào)方法的準(zhǔn)確性；（2）氣象部門(mén)將6小時(shí)降雨量分為6等：0.1—2.5毫米為小雨，2.6—6毫米為中雨，6.1—12毫米為大雨，12.1—25毫米為暴雨，25.1—60毫米為大暴雨，大于60.1毫米為特大暴雨。若按此分級(jí)向公眾預(yù)報(bào)，如何在評(píng)價(jià)方法中考慮公眾的感受？（注：本題數(shù)據(jù)位于壓縮文件C2005Data.rar中,可從/mcm05/problems2005c.asp下載）

10/6/20233例5.1雨量預(yù)報(bào)方法的評(píng)價(jià)（1）請(qǐng)建立數(shù)學(xué)模型來(lái)評(píng)價(jià)兩種6例1雨量預(yù)報(bào)方法的評(píng)價(jià)模型的分析：本題的關(guān)鍵主要是采用Matlab軟件對(duì)所提供的數(shù)據(jù)進(jìn)行分析并得到主要結(jié)論。對(duì)于問(wèn)題一，采用load命令和循環(huán)結(jié)構(gòu)實(shí)現(xiàn)數(shù)據(jù)文件的載入，再?gòu)膶?shí)測(cè)數(shù)據(jù)文件中提出實(shí)測(cè)點(diǎn)位置數(shù)據(jù)，并且依次從預(yù)測(cè)數(shù)據(jù)文件中通過(guò)曲面擬合命令griddata得到相應(yīng)日期、時(shí)段、方法下的對(duì)應(yīng)位置上的預(yù)測(cè)估計(jì)值。然后分別計(jì)算兩種預(yù)報(bào)方法下的實(shí)測(cè)值和預(yù)測(cè)值的偏差并且求對(duì)應(yīng)的總偏差平方和，根據(jù)兩個(gè)總偏差平方和的大小來(lái)得出兩種方法的優(yōu)劣比較，通過(guò)Matlab程序的實(shí)現(xiàn)得到結(jié)果為：第一種方法比第二種方法好。10/6/20234例1雨量預(yù)報(bào)方法的評(píng)價(jià)模型的分析：10/5/20234例1雨量預(yù)報(bào)方法的評(píng)價(jià)對(duì)于問(wèn)題二，我們認(rèn)為公眾的感受主要體現(xiàn)在預(yù)測(cè)與實(shí)際之間偏差的大小程度，然而，實(shí)際測(cè)量值未知，因此我們考慮在分級(jí)預(yù)報(bào)中加入準(zhǔn)確概率的方式來(lái)實(shí)現(xiàn)公眾滿(mǎn)意度的提高。我們認(rèn)為雨量的實(shí)測(cè)值應(yīng)該在預(yù)測(cè)值點(diǎn)處服從正態(tài)分布，通過(guò)合理的假設(shè)和推導(dǎo)，我們得到：正態(tài)分布的均值可以取為雨量預(yù)測(cè)值，不同雨量分級(jí)區(qū)間上的方差可以近似取為對(duì)應(yīng)區(qū)間上的總偏差平方和的平均值。運(yùn)用Matlab軟件編程，可以實(shí)現(xiàn)對(duì)每一個(gè)預(yù)報(bào)數(shù)據(jù)的預(yù)報(bào)內(nèi)容的改變，并且通過(guò)幾個(gè)不同的數(shù)據(jù)體現(xiàn)程序運(yùn)行所得到的結(jié)果。最后，對(duì)模型的缺點(diǎn)進(jìn)行了討論和改進(jìn)。

10/6/20235例1雨量預(yù)報(bào)方法的評(píng)價(jià)對(duì)于問(wèn)題二，我們認(rèn)為公眾的感受主要體例1雨量預(yù)報(bào)方法的評(píng)價(jià)模型的假設(shè)假設(shè)測(cè)量雨量的工具正常不受任何因素的影響，且所得數(shù)據(jù)真實(shí)可靠；測(cè)量點(diǎn)所在位置的等距網(wǎng)格點(diǎn)均視為質(zhì)點(diǎn)；符號(hào)說(shuō)明x(i,j)表示第j種方法下的第i個(gè)實(shí)測(cè)點(diǎn)處的雨量實(shí)測(cè)值y(i,j)表示第j種方法下的第i個(gè)實(shí)測(cè)點(diǎn)處的雨量預(yù)測(cè)估計(jì)值i=1,2,3…….,91*164,j=1,2e(i,j)表示在第j種方法下的第i個(gè)實(shí)測(cè)點(diǎn)處的偏差Y表示實(shí)際值的隨機(jī)變量y*表示預(yù)報(bào)值ST2(j)表示在第j種方法下的所有實(shí)測(cè)點(diǎn)處的總偏差平方和10/6/20236例1雨量預(yù)報(bào)方法的評(píng)價(jià)模型的假設(shè)10/5/20236例1雨量預(yù)報(bào)方法的評(píng)價(jià)問(wèn)題一模型分析：考察下載的lon.dat、lat.dat、020<日期>.SIX及f<日期><i>_dis<j>等數(shù)據(jù)文件。我們可以看到雨量預(yù)報(bào)的網(wǎng)格點(diǎn)有53*47個(gè)。而雨量實(shí)測(cè)點(diǎn)只有91個(gè)，且分布不均勻，因此我們不能直接引用數(shù)據(jù)文件進(jìn)行實(shí)測(cè)數(shù)據(jù)與預(yù)測(cè)數(shù)據(jù)之間偏差的計(jì)算。顯然首先要對(duì)所提供的數(shù)據(jù)文件進(jìn)行處理，由于所引用的文件比較多，而且具有一定的規(guī)律性，所以數(shù)據(jù)文件的載入可用MATLAB命令load和循環(huán)結(jié)構(gòu)來(lái)實(shí)現(xiàn)。當(dāng)有了具體的實(shí)測(cè)點(diǎn)位置對(duì)應(yīng)的實(shí)測(cè)數(shù)據(jù)yi及預(yù)測(cè)數(shù)據(jù)yi*之后,我們可以求出每種方法下所有實(shí)測(cè)點(diǎn)位置上的總偏差平方和。然后比較兩個(gè)總偏差平方和的大小，就可以判斷兩種方法的優(yōu)劣性。

10/6/20237例1雨量預(yù)報(bào)方法的評(píng)價(jià)問(wèn)題一模型分析：10/5/20237例1雨量預(yù)報(bào)方法的評(píng)價(jià)問(wèn)題一模型建立：

通過(guò)Matlab程序編程計(jì)算，具體過(guò)程如下：（1）使用MATLAB命令load和循環(huán)結(jié)構(gòu)將lot.dat、lat.dat、f6181_dis1…、f7304_dis2、020618.six…、020730.six等文件輸入MATLAB程序中備用，其中020618.six…、020730.six等文件系統(tǒng)自動(dòng)在前面添上X。（2）從X020618文件中取第二、第三列作為實(shí)測(cè)點(diǎn)位置。其中第一列對(duì)應(yīng)實(shí)測(cè)點(diǎn)的緯度，第二列對(duì)應(yīng)實(shí)測(cè)點(diǎn)的經(jīng)度。分別用lat、lon中的數(shù)據(jù)作為網(wǎng)格點(diǎn)的橫坐標(biāo)和縱坐標(biāo)。依次取f6181_dis1……,f7304_dis1,f6181_dis2,……,f7304_dis2文件對(duì)應(yīng)點(diǎn)數(shù)據(jù)作為豎坐標(biāo)，張成一個(gè)預(yù)報(bào)值數(shù)據(jù)曲面，用曲面擬合命令griddata得出對(duì)應(yīng)實(shí)測(cè)點(diǎn)處的預(yù)報(bào)估計(jì)值，得到一個(gè)10/6/20238例1雨量預(yù)報(bào)方法的評(píng)價(jià)問(wèn)題一模型建立：10/5/2023例1雨量預(yù)報(bào)方法的評(píng)價(jià)91*328維的對(duì)應(yīng)實(shí)測(cè)點(diǎn)處的預(yù)報(bào)估計(jì)值矩陣yczjz。（3）依次從X020618，…….X020730文件中取第四到第七列組成所有日期，所有時(shí)段下的實(shí)測(cè)數(shù)據(jù)矩陣sczjz。（4）將sczjz分別與yczjz前164列、后164列對(duì)應(yīng)元素相減，得出第一種方法下和第二種方法下所有實(shí)測(cè)點(diǎn)處的偏差，組成矩陣f1pcz和f2pcz。（5）分別求矩陣f1pcz和f2pcz所有元素的平方之和，賦予變量pcpfh(1)和pcpfh(2)。（6）比較變量pcpfh(1)、pcpfh(2)的大小，得出兩種方法優(yōu)劣的比較。各過(guò)程的MATLAB程序?qū)崿F(xiàn)，可分別由M-文件sjzr.m、qsj.m、pcbj.m運(yùn)行得到。

10/6/20239例1雨量預(yù)報(bào)方法的評(píng)價(jià)91*328維的對(duì)應(yīng)實(shí)測(cè)點(diǎn)處的預(yù)報(bào)估例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序sjzr.m：

rqjz=[618:628,701:730];%產(chǎn)生1*41維日期行向量loadlat.DAT;%載入53*47維緯度矩陣loadlon.DAT;%載入53*47維經(jīng)度矩陣forrq=rqjz%用日期行向量作為循環(huán)向量load(['020',int2str(rq),'.SIX'])%載入53*47維實(shí)測(cè)值矩陣,系統(tǒng)自動(dòng)在文件名前加大寫(xiě)字母Xforsd=1:4%指定時(shí)段循環(huán)向量forff=1:2%指定方法循環(huán)向量load(['f',int2str(rq),int2str(sd),'_dis',int2str(ff)])%載入53*47維預(yù)測(cè)值矩陣，共328個(gè)endendend

10/6/202310例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序sjzr.m：10/5/202例1雨量預(yù)報(bào)方法的評(píng)價(jià)將文件全部拷貝在默認(rèn)的work文件加下，運(yùn)行程序sjzr.m：

10/6/202311例1雨量預(yù)報(bào)方法的評(píng)價(jià)將文件全部拷貝在默認(rèn)的work文件加例1雨量預(yù)報(bào)方法的評(píng)價(jià)10/6/202312例1雨量預(yù)報(bào)方法的評(píng)價(jià)10/5/202312例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序huatu.m%畫(huà)一個(gè)時(shí)點(diǎn)預(yù)測(cè)曲面圖與實(shí)測(cè)散點(diǎn)圖比較clf;%預(yù)先清除別的圖象Z=eval(['f',int2str(618),int2str(1),'_dis',int2str(1)]);%取出某時(shí)段預(yù)報(bào)值作為豎坐標(biāo)mesh(lat,lon,Z)%張成預(yù)測(cè)曲面axis([27,36,117,125,0,0.4])shadingflatxlabel('緯度');ylabel('經(jīng)度');title('圖象比較');holdonX1=X020618(:,2);Y1=X020618(:,3);Z1=X020618(:,4);%取出實(shí)測(cè)點(diǎn)數(shù)據(jù)plot3(X1,Y1,Z1,'r*')%畫(huà)出實(shí)測(cè)點(diǎn)數(shù)據(jù)散點(diǎn)圖10/6/202313例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序huatu.m%畫(huà)一個(gè)時(shí)例1雨量預(yù)報(bào)方法的評(píng)價(jià)10/6/202314例1雨量預(yù)報(bào)方法的評(píng)價(jià)10/5/202314例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序qsj.mycr=1;%指定預(yù)測(cè)值矩陣列標(biāo)yczjz=zeros(91,328);%預(yù)先指定預(yù)測(cè)值矩陣為全零矩陣forff=1:2%指定方法循環(huán)變量forrq=rqjz%指定日期循環(huán)變量forsd=1:4%指定時(shí)段循環(huán)變量z0=eval(['f',int2str(rq),int2str(sd),'_dis',int2str(ff)]);%依次取出雨量預(yù)測(cè)值數(shù)據(jù)yczjz(:,ycr)=griddata(lat,lon,z0,X020618(:,2),X020618(:,3));%在張成的預(yù)測(cè)數(shù)據(jù)曲面上擬合對(duì)應(yīng)實(shí)測(cè)點(diǎn)的雨量預(yù)測(cè)值，依次放入預(yù)測(cè)值矩陣對(duì)應(yīng)列ycr=ycr+1;%預(yù)測(cè)值矩陣列標(biāo)向后一列end10/6/202315例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序qsj.m10/5/202315例1雨量預(yù)報(bào)方法的評(píng)價(jià)endendscr=1;%指定實(shí)測(cè)值矩陣列標(biāo)sczjz=zeros(91,164);%預(yù)先指定實(shí)測(cè)值矩陣為全零矩陣forrq=rqjz%指定日期循環(huán)變量sczjz(:,scr:scr+3)=eval(['X020',int2str(rq),'(:,4:7)']);%依次取出雨量實(shí)測(cè)值數(shù)據(jù)放入實(shí)測(cè)值矩陣scr=scr+4;%實(shí)測(cè)值矩陣列標(biāo)向后4列end

10/6/202316例1雨量預(yù)報(bào)方法的評(píng)價(jià)end10/5/2023例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序pcbj.mff1=yczjz(:,1:164);%取出預(yù)測(cè)值矩陣的前164列ff2=yczjz(:,165:328);%取出預(yù)測(cè)值矩陣的后164列f1pcz=sczjz-ff1;%計(jì)算方法1下的偏差f2pcz=sczjz-ff2;%計(jì)算方法2下的偏差pcpfh=[0,0];%預(yù)先產(chǎn)生1*2維偏差平方和零矩陣forpch=1:91%指定偏差行循環(huán)變量forpcr=1:164%指定偏差列循環(huán)變量pcpfh(1)=pcpfh(1)+f1pcz(pch,pcr).^2;%計(jì)算方法1偏差平方和pcpfh(2)=pcpfh(2)+f2pcz(pch,pcr).^2;%計(jì)算方法2偏差平方和10/6/202317例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序pcbj.m10/5/20231例1雨量預(yù)報(bào)方法的評(píng)價(jià)endendpcpfhifpcpfh(1)==pcpfh(2)%比較兩種方法下的偏差平方和display('一樣好')elseifpcpfh(1)<pcpfh(2)display('第一種方法好')elsedisplay('第二種方法好')end首先將數(shù)據(jù)文件減壓縮在matlab程序的work文件夾下，然后依次運(yùn)行M-文件sjzr.m、qsj.m、pcbj.m，得到如下結(jié)論：偏差平方和：226730243880第一種方法比第二種方法好。10/6/202318例1雨量預(yù)報(bào)方法的評(píng)價(jià)end10/5/202318例1雨量預(yù)報(bào)方法的評(píng)價(jià)問(wèn)題二模型分析：在對(duì)問(wèn)題二的分析中，我們只對(duì)方法一相應(yīng)的數(shù)據(jù)進(jìn)行分析。在現(xiàn)實(shí)生活中，公眾對(duì)天氣預(yù)報(bào)的感受主要來(lái)自于對(duì)雨量預(yù)報(bào)等級(jí)和雨量實(shí)際等級(jí)之間的偏差所帶來(lái)的不滿(mǎn)，而不是對(duì)雨量等級(jí)本身的不滿(mǎn)，因此我們認(rèn)為，考慮公眾的感受就是要盡量減少所預(yù)報(bào)等級(jí)和實(shí)際等級(jí)之間偏差的程度。但在預(yù)報(bào)過(guò)程中，我們不可能引用來(lái)自于還未知的實(shí)際數(shù)據(jù)，我們只能從評(píng)價(jià)方法中提供的預(yù)報(bào)數(shù)據(jù)及以前的歷史數(shù)據(jù)來(lái)提高預(yù)報(bào)的準(zhǔn)確性及公眾的感受度。為此我們考慮在分級(jí)預(yù)報(bào)中引入概率。比如：預(yù)報(bào)“中雨”時(shí)，我們可以采取預(yù)報(bào)方式為“中雨，概率為：55%”等。10/6/202319例1雨量預(yù)報(bào)方法的評(píng)價(jià)問(wèn)題二模型分析：10/5/20231例1雨量預(yù)報(bào)方法的評(píng)價(jià)一般的，雨量實(shí)際值是在預(yù)報(bào)值附近振動(dòng)的，而且越接近預(yù)報(bào)值的實(shí)際值的概率越高，離預(yù)報(bào)值越遠(yuǎn)的實(shí)際值的概率越低，因此我們可以認(rèn)為：實(shí)際值服從以預(yù)報(bào)值為中心的正態(tài)分布。假設(shè)實(shí)際值隨機(jī)變量為Y，預(yù)報(bào)值y*則Y~N(y*,σ2)，其中σ2待定?，F(xiàn)假設(shè)對(duì)實(shí)際值Y進(jìn)行n次檢測(cè)得到樣本為Yi，I=1,2,…..,n,對(duì)應(yīng)樣本值為yi，顯然，由概率統(tǒng)計(jì)知識(shí)可知Yi~N(y*,σ2),且是y*的無(wú)偏估計(jì)，是σ2的無(wú)偏估計(jì)。即可用代替y*，s2代替σ2。10/6/202320例1雨量預(yù)報(bào)方法的評(píng)價(jià)一般的，雨量實(shí)際值是在預(yù)報(bào)值附例1雨量預(yù)報(bào)方法的評(píng)價(jià)則：其中ei為每次預(yù)報(bào)的偏差。在此，我們顯然不能對(duì)同一位置進(jìn)行n次偏差檢測(cè)。但是可以用在同一分級(jí)下的實(shí)測(cè)值與預(yù)測(cè)值之間偏差近似代替ei并求σ2值。然后再用各分級(jí)下的σ2值求給定預(yù)報(bào)值所在對(duì)應(yīng)分級(jí)下的預(yù)報(bào)準(zhǔn)確概率P。

10/6/202321例1雨量預(yù)報(bào)方法的評(píng)價(jià)則：10/5/202321例1雨量預(yù)報(bào)方法的評(píng)價(jià)問(wèn)題二模型建立：

由上述分析可知，若將無(wú)雨、小雨、中雨、大雨、暴雨、大暴雨、特大暴雨等7個(gè)等級(jí)所對(duì)應(yīng)的降雨量區(qū)間記為Ij=[aj,bj],則當(dāng)y*∈Ij時(shí)，在第j分級(jí)下的均方差為σ2j，Y~N(y*,σ2j),而且Y的實(shí)際值y∈Ij的概率為：P=預(yù)報(bào)方式為：等級(jí)，概率為：P10/6/202322例1雨量預(yù)報(bào)方法的評(píng)價(jià)問(wèn)題二模型建立：10/5/2023例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序pcbj.mff1=yczjz(:,1:164);%取出預(yù)測(cè)值矩陣的前164列ff2=yczjz(:,165:328);%取出預(yù)測(cè)值矩陣的后164列f1pcz=sczjz-ff1;%計(jì)算方法1下的偏差f2pcz=sczjz-ff2;%計(jì)算方法2下的偏差pcpfh=[0,0];%預(yù)先產(chǎn)生1*2維偏差平方和零矩陣forpch=1:91%指定偏差行循環(huán)變量forpcr=1:164%指定偏差列循環(huán)變量pcpfh(1)=pcpfh(1)+f1pcz(pch,pcr).^2;%計(jì)算方法1偏差平方和pcpfh(2)=pcpfh(2)+f2pcz(pch,pcr).^2;%計(jì)算方法2偏差平方和10/6/202323例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序pcbj.m10/5/20232例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序執(zhí)行過(guò)程描述如下：（1）用if-else-end語(yǔ)句和for循環(huán)對(duì)問(wèn)題一中所得的偏差矩陣進(jìn)行分級(jí)求偏差平方和，并累計(jì)個(gè)數(shù)。（2）計(jì)算各分級(jí)下的偏差平方和與長(zhǎng)度減一的比值，分別作為對(duì)應(yīng)等級(jí)下的正態(tài)分布的均方差的近似值。（3）輸入預(yù)測(cè)值，利用求正態(tài)分布分布函數(shù)命令normcdf及上步中σ值計(jì)算對(duì)應(yīng)分級(jí)中的預(yù)報(bào)準(zhǔn)確概率。相應(yīng)程序見(jiàn)M-文件flpcjz.m、yb.m.10/6/202324例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序執(zhí)行過(guò)程描述如下：10/5/20例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序flpcjz.mwypcf=0;%給‘無(wú)雨’偏差平方和賦初值xypcf=0;%給‘小雨’偏差平方和賦初值z(mì)ypcf=0;%給‘中雨’偏差平方和賦初值dypcf=0;%給‘大雨’偏差平方和賦初值bypcf=0;%給‘暴雨’偏差平方和賦初值dbypcf=0;%給‘大暴雨’偏差平方和賦初值tdypcf=0;%給‘特大暴雨’偏差平方和賦初值wycd=0;%給‘無(wú)雨’偏差個(gè)數(shù)賦初值xycd=0;%給‘小雨’偏差個(gè)數(shù)賦初值z(mì)ycd=0;%給‘中雨’偏差個(gè)數(shù)賦初值dycd=0;%給‘大雨’偏差個(gè)數(shù)賦初值bycd=0;%給‘暴雨’偏差個(gè)數(shù)賦初值10/6/202325例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序flpcjz.m10/5/202例1雨量預(yù)報(bào)方法的評(píng)價(jià)dbycd=0;%給‘大暴雨’偏差個(gè)數(shù)賦初值tdycd=0;%給‘特大暴雨’偏差個(gè)數(shù)賦初值forcol=1:91%用預(yù)測(cè)值矩陣行數(shù)作為循環(huán)變量forrow=1:164,%用預(yù)測(cè)值矩陣列數(shù)作為循環(huán)變量ifyczjz(col,row)<=0.1%判斷預(yù)測(cè)值矩陣對(duì)應(yīng)位置處值的分類(lèi)wypcf=wypcf+f1pcz(col,row).^2;wycd=wycd+1;%按分類(lèi)累加偏差，并且對(duì)應(yīng)個(gè)數(shù)加1elseifyczjz(col,row)>0.1&yczjz(col,row)<=2.5,xypcf=xypcf+f1pcz(col,row).^2;xycd=xycd+1;elseifyczjz(col,row)>2.5&yczjz(col,row)<=6,zypcf=zypcf+f1pcz(col,row).^2;zycd=zycd+1;10/6/202326例1雨量預(yù)報(bào)方法的評(píng)價(jià)dbycd=0;例1雨量預(yù)報(bào)方法的評(píng)價(jià)elseifyczjz(col,row)>6&yczjz(col,row)<=12,dypcf=dypcf+f1pcz(col,row).^2;dycd=dycd+1;elseifyczjz(col,row)>12&yczjz(col,row)<25,bypcf=bypcf+f1pcz(col,row).^2;bycd=bycd+1;elseifyczjz(col,row)>25&yczjz(col,row)<=60,dbypcf=dbypcf+f1pcz(col,row).^2;dbycd=dbycd+1;elsetdypcf=tdypcf+f1pcz(col,row).^2;tdycd=tdycd+1;endendend10/6/202327例1雨量預(yù)報(bào)方法的評(píng)價(jià)elseifyczjz(col,例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序yb.mwz1='無(wú)雨';wz2='小雨';wz3='中雨';wz4='大雨';wz5='暴雨';wz6='大暴雨';wz7='特大暴雨';ifycz<=0.1gl=normcdf(0.1,ycz,sgm(1));wz=wz1;%根據(jù)預(yù)報(bào)值計(jì)算分級(jí)和相應(yīng)概率elseifycz>0.1&ycz<=2.5wz=wz2;10/6/202328例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序yb.m10/5/202328例1雨量預(yù)報(bào)方法的評(píng)價(jià)gl=normcdf(2.5,ycz,sgm(2))-normcdf(0.1,ycz,sgm(2));elseifycz>2.5&ycz<=6gl=normcdf(6,ycz,sgm(3))-normcdf(2.5,ycz,sgm(3));wz=wz3;elseifycz>6&ycz<=12gl=normcdf(12,ycz,sgm(4))-normcdf(6,ycz,sgm(4));wz=wz4;elseifycz>12&ycz<=25gl=normcdf(25,ycz,sgm(5))-normcdf(12,ycz,sgm(5));wz=wz5;elseifycz>25&ycz<=60gl=normcdf(60,ycz,sgm(6))-normcdf(25,ycz,sgm(6));wz=wz6;10/6/202329例1雨量預(yù)報(bào)方法的評(píng)價(jià)gl=normcdf(例1雨量預(yù)報(bào)方法的評(píng)價(jià)elsegl=1-normcdf(60,ycz,sgm(7));wz=wz7;enddisplay([eval('wz'),',概率為:']),display(gl)%輸出等級(jí)和概率

首先運(yùn)行程序flpcjz.m，得到：sgm=[0.08070.60641.95762.74535.191513.5697105.4409]然后：輸入預(yù)報(bào)值：ycz=2.2結(jié)果為：小雨，概率為：0.6893ycz=0.2結(jié)果為：小雨，概率為：0.5654ycz=4.5結(jié)果為：中雨，概率為：0.6248ycz=15.6結(jié)果為：暴雨，概率為：0.720910/6/202330例1雨量預(yù)報(bào)方法的評(píng)價(jià)else10/5/202330例1雨量預(yù)報(bào)方法的評(píng)價(jià)模型的進(jìn)一步討論：

(1)在對(duì)matlab程序進(jìn)行調(diào)試時(shí)，我們觀(guān)察發(fā)現(xiàn)矩陣中出現(xiàn)負(fù)值，但顯然對(duì)天氣預(yù)報(bào)數(shù)據(jù)來(lái)說(shuō)，降雨量是不可能存在負(fù)值的，這個(gè)缺陷是由插值方法的不足所引起的，在程序中我們對(duì)此做了一定的改進(jìn)。見(jiàn)程序qsjxg.m，通過(guò)設(shè)置一個(gè)條件結(jié)構(gòu)，當(dāng)預(yù)測(cè)雨量插值為正時(shí)，保留原值不變；當(dāng)預(yù)測(cè)雨量插值為負(fù)時(shí)，我們?nèi)☆A(yù)測(cè)雨量插值為0。經(jīng)過(guò)此改進(jìn)，我們得到相應(yīng)的問(wèn)題一的答案為：偏差平方和：226690243840第一種方法比第二種方法好。與原答案（偏差平方和：226730243880第一種方法比第二種方法好）相比較，可見(jiàn)插值誤差對(duì)結(jié)果影響不大。10/6/202331例1雨量預(yù)報(bào)方法的評(píng)價(jià)模型的進(jìn)一步討論：10/5/202例1雨量預(yù)報(bào)方法的評(píng)價(jià)(2)觀(guān)察問(wèn)題二的演示答案，我們發(fā)現(xiàn)有些預(yù)報(bào)的概率有接近于50%的情況存在，表明有另外一種雨量等級(jí)概率比較大的可能。為了增加預(yù)報(bào)等級(jí)的可信度，不妨可以采取如下預(yù)報(bào)形式：大到暴雨，大雨概率：55%，暴雨概率40%。只是計(jì)算量有所增加，程序復(fù)雜程度有所增加。10/6/202332例1雨量預(yù)報(bào)方法的評(píng)價(jià)(2)觀(guān)察問(wèn)題二的演示答案，我例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序qsjxg.mycr=1;%指定預(yù)測(cè)值矩陣列標(biāo)yczjz=zeros(91,328);%指定預(yù)測(cè)值矩陣為全零矩陣forff=1:2%指定方法循環(huán)變量forrq=rqjz%指定日期循環(huán)變量forsd=1:4%指定時(shí)段循環(huán)變量z0=eval(['f',int2str(rq),int2str(sd),'_dis',int2str(ff)]);%依次取出雨量預(yù)測(cè)值數(shù)據(jù)yczjz(:,ycr)=griddata(lat,lon,z0,X020618(:,2),X020618(:,3));%在張成的預(yù)測(cè)數(shù)據(jù)曲面上擬合對(duì)應(yīng)實(shí)測(cè)點(diǎn)的雨量預(yù)測(cè)值，依次放入預(yù)測(cè)值矩陣對(duì)應(yīng)列forhb=1:91%指定行標(biāo)循環(huán)變量ifyczjz(hb,ycr)>=0%判斷預(yù)測(cè)擬合值符合條件10/6/202333例1雨量預(yù)報(bào)方法的評(píng)價(jià)程序qsjxg.m10/5/2023例1雨量預(yù)報(bào)方法的評(píng)價(jià)yczjz(hb,ycr)=yczjz(hb,ycr);%符合條件，則值不變elseyczjz(bjl,ycr)=0;%不符合條件，則值為0endendycr=ycr+1;%預(yù)測(cè)值矩陣列標(biāo)向后一列endendendscr=1;%指定實(shí)測(cè)值矩陣列標(biāo)forrq=rqjz%指定日期循環(huán)變量sczjz(:,scr:scr+3)=eval(['X020',int2str(rq),'(:,4:7)']);%依次取出雨量實(shí)測(cè)值數(shù)據(jù)scr=scr+4;%實(shí)測(cè)值矩陣列標(biāo)加4end10/6/202334例1雨量預(yù)報(bào)方法的評(píng)價(jià)yczjz(hb,ycr)=yczj（五）MATLAB應(yīng)用舉例例22004全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽題目飲酒駕車(chē)的數(shù)學(xué)模型飲酒駕車(chē)問(wèn)題主要是分析駕駛員在喝過(guò)一定量的酒后，酒精在體內(nèi)被吸收后，血液中酒精含量上升，影響司機(jī)駕車(chē)，所以司機(jī)飲酒后需經(jīng)過(guò)一段時(shí)間后才能安全駕車(chē)，國(guó)家標(biāo)準(zhǔn)新規(guī)定，車(chē)輛駕駛?cè)藛T血液中的酒精含量大于或等于20毫克/百毫升，小于80毫克/百毫升為飲酒駕車(chē)，血液中酒精含量大于或等于80毫克/百毫升為醉酒駕車(chē)，司機(jī)大李在中午12點(diǎn)喝下一瓶啤酒，6小時(shí)后檢查符合新標(biāo)準(zhǔn)，晚飯地其又喝了一瓶啤酒，他到凌晨2點(diǎn)駕車(chē)，被檢查時(shí)定為飲酒駕車(chē)，為什么喝相同量的酒，兩次結(jié)果不一樣？ 10/6/202335（五）MATLAB應(yīng)用舉例例22004全國(guó)大學(xué)生數(shù)例2飲酒駕車(chē)的數(shù)學(xué)模型討論問(wèn)題：1、對(duì)大李碰到的情況做出合理解釋?zhuān)?、在喝三瓶啤酒或半斤白酒后多長(zhǎng)時(shí)間內(nèi)駕車(chē)會(huì)違反標(biāo)準(zhǔn)，喝酒時(shí)間長(zhǎng)短不同情況會(huì)怎樣？3、分析當(dāng)司機(jī)喝酒后何時(shí)血液中的酒精含量最高；4、如果該司機(jī)想天天喝酒還能否開(kāi)車(chē)；10/6/202336例2飲酒駕車(chē)的數(shù)學(xué)模型討論問(wèn)題：10/5/202336例2飲酒駕車(chē)的數(shù)學(xué)模型摘要本文解決的是一個(gè)司機(jī)安全駕車(chē)與飲酒的問(wèn)題，目的是通過(guò)建立一個(gè)數(shù)學(xué)模型（結(jié)合新的國(guó)家駕駛員飲酒標(biāo)準(zhǔn)）分析司機(jī)如何適量飲酒不會(huì)影響正常的安全駕駛。根據(jù)一定合理的假設(shè)，建立人體內(nèi)酒精濃度隨時(shí)間變化的微分方程模型，并通過(guò)擬合曲線(xiàn)對(duì)數(shù)據(jù)進(jìn)行分析。在不同飲酒方式下進(jìn)行分類(lèi)討論，得出體內(nèi)酒精濃度隨時(shí)間的變化函數(shù)。在討論過(guò)程中，我們得到兩個(gè)結(jié)論：在短時(shí)間喝酒形式下，達(dá)到最大值的時(shí)間為1.23小時(shí)，與喝酒量無(wú)關(guān)；在長(zhǎng)時(shí)間喝酒形式下，喝酒結(jié)束時(shí)酒精含量最高。10/6/202337例2飲酒駕車(chē)的數(shù)學(xué)模型摘要10/5/202337例2飲酒駕車(chē)的數(shù)學(xué)模型模型假設(shè)1、酒精從胃轉(zhuǎn)移到體液的速率與胃中的酒精濃度成正比。2、酒精從體液轉(zhuǎn)移到體外的速率與體液中的酒精濃度成正比。3、酒精從胃轉(zhuǎn)移到體液的過(guò)程中沒(méi)有損失。4、測(cè)量設(shè)備完善，不考慮不同因素所造成的誤差。5、酒精在體液中均勻分布。10/6/202338例2飲酒駕車(chē)的數(shù)學(xué)模型模型假設(shè)10/5/202338例2飲酒駕車(chē)的數(shù)學(xué)模型符號(hào)說(shuō)明k:酒精從體外進(jìn)入胃的速率；f1(t):酒精從胃轉(zhuǎn)移到體液的速率；f2(t):酒精從體液轉(zhuǎn)移到體外的速率；X(t):胃里的酒精含量；Y(t)：體液中酒精含量；V0:體液的容積；K1:酒精從胃轉(zhuǎn)移到體液的轉(zhuǎn)移速率系數(shù)；K2:酒精從體液轉(zhuǎn)移到體外的轉(zhuǎn)移速率系數(shù)；C(t):體液中的酒精濃度。10/6/202339例2飲酒駕車(chē)的數(shù)學(xué)模型符號(hào)說(shuō)明10/5/202339例2飲酒駕車(chē)的數(shù)學(xué)模型D0：短時(shí)間喝酒情況下進(jìn)入胃中的初始酒精量。T：較長(zhǎng)時(shí)間喝酒所用的時(shí)間或達(dá)到濃度最大值所需時(shí)間。模型分析：假設(shè)酒精先以速率k0進(jìn)入胃中，然后以速率f1(t)從胃進(jìn)入體液，再以速率f2(t)從體液中排到體外。根據(jù)假設(shè)可以建立如圖所示的帶有吸收室的單房室系統(tǒng)，其中胃為吸收室，體液為中心室。10/6/202340例2飲酒駕車(chē)的數(shù)學(xué)模型D0：短時(shí)間喝酒情況下進(jìn)入胃中的初始酒例2飲酒駕車(chē)的數(shù)學(xué)模型根據(jù)酒精從胃進(jìn)入體液的速度f(wàn)1(t)與胃中的酒精量成正比，速率系數(shù)為K1；酒精從血液中排出的速率f2(t)與血液中的酒精量y(t)成正比，速率系數(shù)為K2，可以建立方程如下：10/6/202341例2飲酒駕車(chē)的數(shù)學(xué)模型根據(jù)酒精從胃進(jìn)入體液的速度f(wàn)1(t)與例2飲酒駕車(chē)的數(shù)學(xué)模型10/6/202342例2飲酒駕車(chē)的數(shù)學(xué)模型10/5/202342例2飲酒駕車(chē)的數(shù)學(xué)模型利用MATLAB的dsolve()函數(shù)，編輯：clear;symsxyy0k0k1k2;[x,y]=dsolve('Dx=k0-k1*x,Dy=k1*x-k2*y',…'x(0)=x0,y(0)=y0','t')pretty(simple(x)),pretty(simple(y))得到：10/6/202343例2飲酒駕車(chē)的數(shù)學(xué)模型利用MATLAB的dsolve()函數(shù)例2飲酒駕車(chē)的數(shù)學(xué)模型10/6/202344例2飲酒駕車(chē)的數(shù)學(xué)模型10/5/202344例2飲酒駕車(chē)的數(shù)學(xué)模型模型討論：

當(dāng)酒是在較短時(shí)間內(nèi)喝時(shí)此時(shí)有x(0)=x0=D0，k0=0,y(0)=y0=0:又酒精濃度為酒精量與體液容積之比，即：10/6/2