MATLAB使用教程(5)应用举例课件.ppt

1、2022-6-16黄建华制作0（五） MATLAB应用举例v例5.1 2005高教社杯全国大学生数学建模竞赛题目 C雨量雨量预报方法的评价预报方法的评价 雨量预报对农业生产和城市工作和生活有重要作用，但准确、及时地对雨量作出预报是一个十分困难的问题，广受世界各国关注。我国某地气象台和气象研究所正在研究6小时雨量预报方法，即每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度、北纬32度附近的5347的等距网格点上。同时设立91个观测站点实测这些时段的实际雨量，由于各种条件的限制，站点的设置是不均匀的

2、。 2022-6-16黄建华制作1例5.1 雨量预报方法的评价雨量预报方法的评价气象部门希望建立一种科学评价预报方法好坏的数学模型与方法。气象部门提供了41天的用两种不同方法的预报数据和相应的实测数据。预报数据在文件夹FORECAST中，实测数据在文件夹MEASURING中，其中的文件都可以用Windows系统的“写字板”程序打开阅读。FORECAST中的文件lon.dat和lat.dat分别包含网格点的经纬度，其余文件名为_dis1和_dis2，例如f6181_dis1中包含2002年6月18日晚上20点采用第一种方法预报的第一时段数据（其2491个数据为该时段各网格点的雨量），而f6183

3、_dis2中包含2002年6月18日晚上20点采用第二种方法预报的第三时段数据。2022-6-16黄建华制作2例5.1 雨量预报方法的评价雨量预报方法的评价 MEASURING中包含了41个名为.SIX的文件，如020618.SIX表示2002年6月18日晚上21点开始的连续4个时段各站点的实测数据（雨量），这些文件的数据格式是：站号 纬度 经度 第1段 第2段 第3段 第4段 58138 32.9833 118.5167 0.0000 0.2000 10.1000 3.100058139 33.3000 118.8500 0.0000 0.0000 4.6000 7.400058141 33

4、.6667 119.2667 0.0000 0.0000 1.1000 1.400058143 33.8000 119.8000 0.0000 0.0000 0.0000 1.800058146 33.4833 119.8167 0.0000 0.0000 1.5000 1.9000雨量用毫米做单位，小于0.1毫米视为无雨。 2022-6-16黄建华制作3例5.1 雨量预报方法的评价雨量预报方法的评价（1）请建立数学模型来评价两种6小时雨量预报方法的准确性；（2）气象部门将6小时降雨量分为6等：0.12.5毫米为小雨，2.66毫米为中雨，6.112毫米为大雨，12.125毫米为暴雨，25.16

5、0毫米为大暴雨，大于60.1毫米为特大暴雨。若按此分级向公众预报，如何在评价方法中考虑公众的感受？ （注：本题数据位于压缩文件C2005Data.rar中, 可从http:/ 2022-6-16黄建华制作4例5.1 雨量预报方法的评价雨量预报方法的评价v模型的分析：模型的分析：本题的本题的关键主要是采用关键主要是采用Matlab软件对所提供的数据进软件对所提供的数据进行分析并得到主要结论。行分析并得到主要结论。对于问题一，采用对于问题一，采用load命令和命令和循环结构循环结构实现数据文件实现数据文件的载入，再从实测数据文件中提出实测点位置数据，的载入，再从实测数据文件中提出实测点位置数据，并

6、且依次从预测数据文件中通过曲面拟合命令并且依次从预测数据文件中通过曲面拟合命令griddata得到相应日期、时段、方法下的对应位置上得到相应日期、时段、方法下的对应位置上的预测估计值。然后分别计算两种预报方法下的实测的预测估计值。然后分别计算两种预报方法下的实测值和预测值的偏差并且求对应的总偏差平方和，根据值和预测值的偏差并且求对应的总偏差平方和，根据两个总偏差平方和的大小来得出两种方法的优劣比较，两个总偏差平方和的大小来得出两种方法的优劣比较，通过通过Matlab程序的实现得到结果为：第一种方法比第程序的实现得到结果为：第一种方法比第二种方法好。二种方法好。2022-6-16黄建华制作5例5

7、.1 雨量预报方法的评价雨量预报方法的评价对于问题二，我们认为公众的感受主要体现在预测与对于问题二，我们认为公众的感受主要体现在预测与实际之间偏差的大小程度，然而，实际测量值未知，实际之间偏差的大小程度，然而，实际测量值未知，因此我们考虑在分级预报中加入准确概率的方式来实因此我们考虑在分级预报中加入准确概率的方式来实现公众满意度的提高。我们认为雨量的实测值应该在现公众满意度的提高。我们认为雨量的实测值应该在预测值点处服从正态分布，通过合理的假设和推导，预测值点处服从正态分布，通过合理的假设和推导，我们得到：正态分布的均值可以取为雨量预测值，不我们得到：正态分布的均值可以取为雨量预测值，不同雨量

8、分级区间上的方差可以近似取为对应区间上的同雨量分级区间上的方差可以近似取为对应区间上的总偏差平方和的平均值。运用总偏差平方和的平均值。运用Matlab软件编程，可以软件编程，可以实现对每一个预报数据的预报内容的改变，并且通过实现对每一个预报数据的预报内容的改变，并且通过几个不同的数据体现程序运行所得到的结果。最后，几个不同的数据体现程序运行所得到的结果。最后，对模型的缺点进行了讨论和改进。对模型的缺点进行了讨论和改进。 2022-6-16黄建华制作6例5.1 雨量预报方法的评价雨量预报方法的评价v模型的假设模型的假设假设测量雨量的工具正常不受任何因素的影响，且所得数据真假设测量雨量的工具正常不

9、受任何因素的影响，且所得数据真实可靠；实可靠；测量点所在位置的等距网格点均视为质点；测量点所在位置的等距网格点均视为质点；v符号说明符号说明x(i,j) 表示第表示第j种方法下的第种方法下的第i 个实测点处的雨量实测值个实测点处的雨量实测值y(i,j) 表示第表示第j种方法下的第种方法下的第i 个实测点处的雨量预测估计值个实测点处的雨量预测估计值 i=1,2,3.,91*164,j=1,2e(i,j) 表示在第表示在第j种方法下的第种方法下的第i 个实测点处的偏差个实测点处的偏差Y 表示实际值的随机变量表示实际值的随机变量y* 表示预报值表示预报值ST2(j) 表示在第表示在第j种方法下的所有

10、实测点处的总偏差平方和种方法下的所有实测点处的总偏差平方和2022-6-16黄建华制作7例5.1 雨量预报方法的评价雨量预报方法的评价v问题一模型分析：问题一模型分析： 考察下载的考察下载的lon.dat 、 lat.dat、020.SIX及及f_dis等数据文件。我们可以看到雨量预报的网格点有等数据文件。我们可以看到雨量预报的网格点有53*47个个。而雨量实测点只有。而雨量实测点只有91个个，且，且分布不均匀分布不均匀，因此我们，因此我们不能直接引用数据文件进行实测数据与预测数据之间偏差的计不能直接引用数据文件进行实测数据与预测数据之间偏差的计算。显然首先要对所提供的数据文件进行处理，由于所

11、引用的算。显然首先要对所提供的数据文件进行处理，由于所引用的文件比较多，而且具有一定的规律性，所以数据文件的载入可文件比较多，而且具有一定的规律性，所以数据文件的载入可用用MATLAB命令命令 load和循环结构来实现。当有了具体的实测和循环结构来实现。当有了具体的实测点位置对应的实测数据点位置对应的实测数据yi及预测数据及预测数据yi*之后之后,我们可以求出每我们可以求出每种方法下所有实测点位置种方法下所有实测点位置上的总偏差平方和。然后比较两个总上的总偏差平方和。然后比较两个总偏差平方和的大小，就可以判断两种方法的优劣性。偏差平方和的大小，就可以判断两种方法的优劣性。 2022-6-16黄

12、建华制作8例5.1 雨量预报方法的评价雨量预报方法的评价v问题一模型建立： 通过通过Matlab程序编程计算，具体过程如下：程序编程计算，具体过程如下：（1）使用）使用MATLAB命令命令 load和循环结构将和循环结构将lot.dat、lat.dat、f6181_dis1、f7304_dis2、020618.six、020730.six等等文件输入文件输入MATLAB程序中备用，其中程序中备用，其中020618.six、020730.six等文件系统自动在前面添上等文件系统自动在前面添上X。（2）从）从X020618文件中取第二、第三列作为实测点位置。其文件中取第二、第三列作为实测点位置。其

13、中第一列对应实测点的纬度，第二列对应实测点的经度。分别中第一列对应实测点的纬度，第二列对应实测点的经度。分别用用lat、lon中的数据作为网格点的横坐标和纵坐标。依次取中的数据作为网格点的横坐标和纵坐标。依次取f6181_dis1,f7304_dis1,f6181_dis2,f7304_dis2文文件对应点数据作为竖坐标，张成一个预报值数据曲面，用曲面件对应点数据作为竖坐标，张成一个预报值数据曲面，用曲面拟合命令拟合命令griddata得出对应实测点处的预报估计值，得到一个得出对应实测点处的预报估计值，得到一个2022-6-16黄建华制作9例5.1 雨量预报方法的评价雨量预报方法的评价91*3

14、28维的对应实测点处的预报估计值矩阵维的对应实测点处的预报估计值矩阵yczjz。（3）依次从）依次从X020618，.X020730文件中取第四到第七列文件中取第四到第七列组成所有日期，所有时段下的实测数据矩阵组成所有日期，所有时段下的实测数据矩阵sczjz。（4）将）将sczjz分别与分别与yczjz前前164列、后列、后164列对应元素相减，列对应元素相减，得出第一种方法下和第二种方法下所有实测点处的偏差，组成得出第一种方法下和第二种方法下所有实测点处的偏差，组成矩阵矩阵f1pcz和和f2pcz。 （5）分别求矩阵）分别求矩阵f1pcz和和f2pcz所有元素的平方之和，赋予变所有元素的平方

15、之和，赋予变量量pcpfh(1)和和 pcpfh(2)。（6）比较变量）比较变量pcpfh(1)、 pcpfh(2)的大小，得出两种方法优的大小，得出两种方法优劣的比较。劣的比较。各过程的各过程的MATLAB程序实现，可分别由程序实现，可分别由M-文件文件sjzr.m、qsj.m、pcbj.m运行得到。运行得到。 2022-6-16黄建华制作10例5.1 雨量预报方法的评价雨量预报方法的评价v程序sjzr.m： rqjz=618:628,701:730; %产生产生1*41维日期行向量维日期行向量load lat.DAT; %载入载入53*47维纬度矩阵维纬度矩阵load lon.DAT; %

16、载入载入53*47维经度矩阵维经度矩阵for rq=rqjz %用日期行向量作为循环向量用日期行向量作为循环向量 load(020,int2str(rq),.SIX) %载入载入53*47维实测值矩阵维实测值矩阵,系系 统自动在文件名前加大写字母统自动在文件名前加大写字母X for sd=1:4 %指定时段循环向量指定时段循环向量 for ff=1:2 %指定方法循环向量指定方法循环向量 load(f,int2str(rq),int2str(sd),_dis,int2str(ff) %载入载入53*47维预测值矩阵，共维预测值矩阵，共328个个 end endend 2022-6-16黄建华制

17、作11例5.1 雨量预报方法的评价雨量预报方法的评价将文件全部拷贝在默认的work文件加下，运行程序sjzr.m： 2022-6-16黄建华制作12例5.1 雨量预报方法的评价雨量预报方法的评价2022-6-16黄建华制作13例5.1 雨量预报方法的评价雨量预报方法的评价v程序huatu.m %画一个时点预测曲面图与实测散点图比较画一个时点预测曲面图与实测散点图比较clf; %预先清除别的图象预先清除别的图象 Z=eval(f,int2str(618),int2str(1),_dis,int2str(1); %取出某时段预报值作为竖坐标取出某时段预报值作为竖坐标mesh(lat,lon,Z)

18、%张成预测曲面张成预测曲面 axis(27,36,117,125,0,0.4)shading flat xlabel(纬度纬度); ylabel(经度经度);title(图象比较图象比较);hold on X1=X020618(:,2);Y1= X020618(:,3); Z1= X020618 (:,4); %取出实测点数据取出实测点数据plot3(X1,Y1,Z1,r*) %画出实测点数据散点图画出实测点数据散点图2022-6-16黄建华制作14例5.1 雨量预报方法的评价雨量预报方法的评价2022-6-16黄建华制作15例5.1 雨量预报方法的评价雨量预报方法的评价v程序qsj.mycr

19、=1; %指定预测值矩阵列标指定预测值矩阵列标yczjz=zeros(91,328); %预先指定预测值矩阵为全零矩阵预先指定预测值矩阵为全零矩阵for ff=1:2 %指定方法循环变量指定方法循环变量 for rq=rqjz %指定日期循环变量指定日期循环变量 for sd=1:4 %指定时段循环变量指定时段循环变量 z0=eval(f,int2str(rq),int2str(sd),_dis,int2str(ff); %依次取出雨量预测值数据依次取出雨量预测值数据 yczjz(:,ycr)=griddata(lat,lon,z0,X020618(:,2),X020618(:,3); %在张

20、成的预测数据曲面上拟合对应实测点的雨量预测在张成的预测数据曲面上拟合对应实测点的雨量预测 值，依次放入预测值矩阵对应列值，依次放入预测值矩阵对应列 ycr=ycr+1; %预测值矩阵列标向后一列预测值矩阵列标向后一列 end 2022-6-16黄建华制作16例5.1 雨量预报方法的评价雨量预报方法的评价 endendscr=1; %指定实测值矩阵列标指定实测值矩阵列标sczjz=zeros(91,164); %预先指定实测值矩阵为全零矩阵预先指定实测值矩阵为全零矩阵for rq=rqjz %指定日期循环变量指定日期循环变量 sczjz(:,scr:scr+3)=eval(X020,int2st

21、r(rq),(:,4:7); %依次取出雨量实测值数据放入实测值矩阵依次取出雨量实测值数据放入实测值矩阵 scr=scr+4; %实测值矩阵列标向后实测值矩阵列标向后4列列 end 2022-6-16黄建华制作17例5.1 雨量预报方法的评价雨量预报方法的评价v程序pcbj.mff1=yczjz(:,1:164); %取出预测值矩阵的前取出预测值矩阵的前164列列ff2=yczjz(:,165:328); %取出预测值矩阵的后取出预测值矩阵的后164列列f1pcz=sczjz-ff1; %计算方法计算方法1下的偏差下的偏差f2pcz=sczjz-ff2; %计算方法计算方法2下的偏差下的偏差p

22、cpfh=0,0; %预先产生预先产生1*2维偏差平方和零矩阵维偏差平方和零矩阵for pch=1:91 %指定偏差行循环变量指定偏差行循环变量 for pcr=1:164 %指定偏差列循环变量指定偏差列循环变量 pcpfh(1)=pcpfh(1)+f1pcz(pch,pcr).2; %计算方法计算方法1偏差平方和偏差平方和 pcpfh(2)=pcpfh(2)+f2pcz(pch,pcr).2; %计算方法计算方法2偏差平方和偏差平方和2022-6-16黄建华制作18例5.1 雨量预报方法的评价雨量预报方法的评价 endendpcpfhif pcpfh(1)=pcpfh(2) %比较两种方法下

23、的偏差平方和比较两种方法下的偏差平方和 display(一样好一样好)else if pcpfh(1)pcpfh(2) display(第一种方法好第一种方法好)else display(第二种方法好第二种方法好)endv首先将数据文件减压缩在首先将数据文件减压缩在matlab程序的程序的work文件夹下，然后依次运行文件夹下，然后依次运行M-文件文件sjzr.m、qsj.m、pcbj.m，得到如下结论：偏差平方和：，得到如下结论：偏差平方和：226730 243880 第一种方法比第二种方法好。第一种方法比第二种方法好。2022-6-16黄建华制作19例5.1 雨量预报方法的评价雨量预报方法

24、的评价v问题二模型分析：问题二模型分析： 在对问题二的分析中，我们只对方法一相应的数据进行分在对问题二的分析中，我们只对方法一相应的数据进行分析。析。在现实生活中，公众对天气预报的感受主要来自于对雨量预报在现实生活中，公众对天气预报的感受主要来自于对雨量预报等级和雨量实际等级之间的偏差所带来的不满，而不是对雨量等级和雨量实际等级之间的偏差所带来的不满，而不是对雨量等级本身的不满，因此我们认为，考虑公众的感受就是要尽量等级本身的不满，因此我们认为，考虑公众的感受就是要尽量减少所预报等级和实际等级之间偏差的程度。但在预报过程中，减少所预报等级和实际等级之间偏差的程度。但在预报过程中，我们不可能引用

25、来自于还未知的实际数据，我们只能从评价方我们不可能引用来自于还未知的实际数据，我们只能从评价方法中提供的预报数据及以前的历史数据来提高预报的准确性及法中提供的预报数据及以前的历史数据来提高预报的准确性及公众的感受度。为此我们考虑在分级预报中引入概率。比如：公众的感受度。为此我们考虑在分级预报中引入概率。比如：预报预报“中雨中雨”时，我们可以采取预报方式为时，我们可以采取预报方式为“中雨，概率为：中雨，概率为：55%”等。等。2022-6-16黄建华制作20例5.1 雨量预报方法的评价雨量预报方法的评价 一般的，雨量实际值是在预报值附近振动的，而且越接近预一般的，雨量实际值是在预报值附近振动的，

26、而且越接近预报值的实际值的概率越高，离预报值越远的实际值的概率越低，报值的实际值的概率越高，离预报值越远的实际值的概率越低，因此我们可以认为：实际值服从以预报值为中心的正态分布。因此我们可以认为：实际值服从以预报值为中心的正态分布。假设实际值随机变量为假设实际值随机变量为Y，预报值，预报值y* 则则YN(y*, 2 2 )，其中，其中 2 2待定。现假设对实际值待定。现假设对实际值Y进行进行n次检测得到样本为次检测得到样本为Yi，I=1,2,.,n,对应样本值为对应样本值为yi，显然，由概率统计知识可知，显然，由概率统计知识可知YiN(y*, 2 2 ),且且 是是y*的无偏估计，的无偏估计，

27、 是是 2 2 的无偏估计。即可用的无偏估计。即可用 代替代替y*， s2代替代替 2 2。_1iyyn2_211()1niisyyn_1iyyn2022-6-16黄建华制作21例5.1 雨量预报方法的评价雨量预报方法的评价则则 ： 其中其中ei为每次预报的偏差。为每次预报的偏差。在此，我们显然不能对同一位置进行在此，我们显然不能对同一位置进行n次偏差检测。但是可以次偏差检测。但是可以用在同一分级下的实测值与预测值之间偏差近似代替用在同一分级下的实测值与预测值之间偏差近似代替ei并求并求2 2值。然后再用各分级下的值。然后再用各分级下的2 2值求给定预报值所在对应分级值求给定预报值所在对应分级

28、下的预报准确概率下的预报准确概率P。 222*21111()11nniiiisyyenn2022-6-16黄建华制作22例5.1 雨量预报方法的评价雨量预报方法的评价v问题二模型建立： 由上述分析可知，若将无雨、小雨、中雨、大雨、暴雨、大暴由上述分析可知，若将无雨、小雨、中雨、大雨、暴雨、大暴雨、特大暴雨等雨、特大暴雨等7个等级所对应的降雨量区间记为个等级所对应的降雨量区间记为Ij=aj,bj, 则当则当y*Ij时，在第时，在第j分级下的均方差为分级下的均方差为2 2 j，YN(y*, 2 2 j),而且而且Y的实际值的实际值yIj的概率为：的概率为：P= 预报方式为：等级，概率为：预报方式为

29、：等级，概率为：P*22()212jjxybaedx2022-6-16黄建华制作23例5.1 雨量预报方法的评价雨量预报方法的评价v程序pcbj.mff1=yczjz(:,1:164); %取出预测值矩阵的前取出预测值矩阵的前164列列ff2=yczjz(:,165:328); %取出预测值矩阵的后取出预测值矩阵的后164列列f1pcz=sczjz-ff1; %计算方法计算方法1下的偏差下的偏差f2pcz=sczjz-ff2; %计算方法计算方法2下的偏差下的偏差pcpfh=0,0; %预先产生预先产生1*2维偏差平方和零矩阵维偏差平方和零矩阵for pch=1:91 %指定偏差行循环变量指定

30、偏差行循环变量 for pcr=1:164 %指定偏差列循环变量指定偏差列循环变量 pcpfh(1)=pcpfh(1)+f1pcz(pch,pcr).2; %计算方法计算方法1偏差平方和偏差平方和 pcpfh(2)=pcpfh(2)+f2pcz(pch,pcr).2; %计算方法计算方法2偏差平方和偏差平方和2022-6-16黄建华制作24例5.1 雨量预报方法的评价雨量预报方法的评价程序执行过程描述如下：程序执行过程描述如下：（1）用）用if-else-end语句和语句和for循环对问题一中所得的偏差矩阵循环对问题一中所得的偏差矩阵进行分级求偏差平方和，并累计个数。进行分级求偏差平方和，并累

31、计个数。（2）计算各分级下的偏差平方和与长度减一的比值，分别作）计算各分级下的偏差平方和与长度减一的比值，分别作为对应等级下的正态分布的均方差为对应等级下的正态分布的均方差 的近似值。的近似值。（3）输入预测值，利用求正态分布分布函数命令）输入预测值，利用求正态分布分布函数命令normcdf及及上步中上步中 值计算对应分级中的预报准确概率。值计算对应分级中的预报准确概率。相应程序见相应程序见M-文件文件flpcjz.m、yb.m.2022-6-16黄建华制作25例5.1 雨量预报方法的评价雨量预报方法的评价v程序flpcjz.mwypcf=0; %给给无雨无雨偏差平方和赋初值偏差平方和赋初值x

32、ypcf=0; %给给小雨小雨偏差平方和赋初值偏差平方和赋初值zypcf=0; %给给中雨中雨偏差平方和赋初值偏差平方和赋初值dypcf=0; %给给大雨大雨偏差平方和赋初值偏差平方和赋初值bypcf=0; %给给暴雨暴雨偏差平方和赋初值偏差平方和赋初值dbypcf=0; %给给大暴雨大暴雨偏差平方和赋初值偏差平方和赋初值tdypcf=0; %给给特大暴雨特大暴雨偏差平方和赋初值偏差平方和赋初值wycd=0; %给给无雨无雨偏差个数赋初值偏差个数赋初值xycd=0; %给给小雨小雨偏差个数赋初值偏差个数赋初值zycd=0; %给给中雨中雨偏差个数赋初值偏差个数赋初值dycd=0; %给给大雨大

33、雨偏差个数赋初值偏差个数赋初值bycd=0; %给给暴雨暴雨偏差个数赋初值偏差个数赋初值2022-6-16黄建华制作26例5.1 雨量预报方法的评价雨量预报方法的评价dbycd=0; %给给大暴雨大暴雨偏差个数赋初值偏差个数赋初值tdycd=0; %给给特大暴雨特大暴雨偏差个数赋初值偏差个数赋初值for col=1:91 %用预测值矩阵行数作为循环变量用预测值矩阵行数作为循环变量 for row=1:164, %用预测值矩阵列数作为循环变量用预测值矩阵列数作为循环变量 if yczjz(col,row)0.1&yczjz(col,row)2.5&yczjz(col,row)6&yczjz(co

34、l,row)12&yczjz(col,row)25&yczjz(col,row)=60, dbypcf=dbypcf +f1pcz(col,row).2; dbycd=dbycd+1; else tdypcf=tdypcf +f1pcz(col,row).2; tdycd=tdycd+1;end end end2022-6-16黄建华制作28例5.1 雨量预报方法的评价雨量预报方法的评价v程序yb.mwz1=无雨无雨 ;wz2=小雨小雨 ;wz3=中雨中雨 ;wz4=大雨大雨 ;wz5=暴雨暴雨 ;wz6=大暴雨大暴雨 ;wz7=特大暴雨特大暴雨;if ycz0.1&ycz2.5&ycz6&y

35、cz12&ycz25&ycz=0 %判断预测拟合值符合条件判断预测拟合值符合条件2022-6-16黄建华制作34例5.1 雨量预报方法的评价雨量预报方法的评价yczjz(hb,ycr)=yczjz(hb,ycr); %符合条件，则值不变符合条件，则值不变else yczjz(bjl,ycr)=0; %不符合条件，则值为不符合条件，则值为0endendycr=ycr+1; %预测值矩阵列标向后一列预测值矩阵列标向后一列endendendscr=1; %指定实测值矩阵列标指定实测值矩阵列标for rq=rqjz %指定日期循环变量指定日期循环变量 sczjz(:,scr:scr+3)=eval(X

36、020,int2str(rq),(:,4:7); %依次取出雨量实测值数据依次取出雨量实测值数据 scr=scr+4;%实测值矩阵列标加实测值矩阵列标加4end2022-6-16黄建华制作35（五） MATLAB应用举例v例5.2 2004全国大学生数学建模竞赛题目 饮酒驾车的数学模型饮酒驾车的数学模型 饮酒驾车问题主要是分析驾驶员在喝过一定量的酒后，酒精在体内被吸收后，血液中酒精含量上升，影响司机驾车，所以司机饮酒后需经过一段时间后才能安全驾车，国家标准新规定，车辆驾驶人员血液中的酒精含量大于或等于20毫克/百毫升，小于80毫克/百毫升为饮酒驾车，血液中酒精含量大于或等于80毫克/百毫升为醉

37、酒驾车，司机大李在中午12点喝下一瓶啤酒，6小时后检查符合新标准，晚饭地其又喝了一瓶啤酒，他到凌晨2点驾车，被检查时定为饮酒驾车，为什么喝相同量的酒，两次结果不一样？2022-6-16黄建华制作36例5.2饮酒驾车的数学模型饮酒驾车的数学模型讨论问题： 1、对大李碰到的情况做出合理解释；2、在喝三瓶啤酒或半斤白酒后多长时间内驾车会违反标准，喝酒时间长短不同情况会怎样？3、分析当司机喝酒后何时血液中的酒精含量最高；4、如果该司机想天天喝酒还能否开车；2022-6-16黄建华制作37例5.2饮酒驾车的数学模型饮酒驾车的数学模型v摘要 本文解决的是一个司机安全驾车与饮酒的问题，目的是通过建立一个数学

38、模型（结合新的国家驾驶员饮酒标准）分析司机如何适量饮酒不会影响正常的安全驾驶。根据一定合理的假设，建立人体内酒精浓度随时间变化的微分方程模型，并通过拟合曲线对数据进行分析。在不同饮酒方式下进行分类讨论，得出体内酒精浓度随时间的变化函数。在讨论过程中，我们得到两个结论：在短时间喝酒形式下，达到最大值的时间为1.23小时，与喝酒量无关；在长时间喝酒形式下，喝酒结束时酒精含量最高。 2022-6-16黄建华制作38例5.2饮酒驾车的数学模型饮酒驾车的数学模型v模型假设1、酒精从胃转移到体液的速率与胃中的酒精浓度成正比。2、酒精从体液转移到体外的速率与体液中的酒精浓度成正比。3、酒精从胃转移到体液的过

39、程中没有损失。4、测量设备完善，不考虑不同因素所造成的误差。5、酒精在体液中均匀分布。2022-6-16黄建华制作39例5.2饮酒驾车的数学模型饮酒驾车的数学模型v符号说明k :酒精从体外进入胃的速率；f 1(t):酒精从胃转移到体液的速率；f 2(t):酒精从体液转移到体外的速率；X(t):胃里的酒精含量；Y(t)：体液中酒精含量；V 0:体液的容积；K 1:酒精从胃转移到体液的转移速率系数；K2 :酒精从体液转移到体外的转移速率系数；C(t):体液中的酒精浓度。2022-6-16黄建华制作40例5.2饮酒驾车的数学模型饮酒驾车的数学模型D0：短时间喝酒情况下进入胃中的初始酒精量。T：较长时

40、间喝酒所用的时间或达到浓度最大值所需时间。v模型分析： 假设酒精先以速率k0进入胃中，然后以速率f1(t)从胃进入体液，再以速率f2(t)从体液中排到体外。根据假设可以建立如图所示的带有吸收室的单房室系统，其中胃为吸收室，体液为中心室。 2022-6-16黄建华制作41例5.2饮酒驾车的数学模型饮酒驾车的数学模型根据酒精从胃进入体液的速度f1(t)与胃中的酒精量成正比，速率系数为K1；酒精从血液中排出的速率f2(t)与血液中的酒精量y(t)成正比，速率系数为K2，可以建立方程如下： 2022-6-16黄建华制作42例5.2饮酒驾车的数学模型饮酒驾车的数学模型)()(11txktf)()(22t

41、yktf)()(10tfkdttdx)()(10txkkdttdx021)0()()()(yytftfdttdy)()()(21tyktxkdttdy2022-6-16黄建华制作43例5.2饮酒驾车的数学模型饮酒驾车的数学模型002110)0(,)0()()()()()(yyxxtyktxkdttdytxkkdttdx利用MATLAB的dsolve()函数，编辑： clear; syms x y y0 k0 k1 k2;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

42、)得到：2022-6-16黄建华制作44例5.2饮酒驾车的数学模型饮酒驾车的数学模型tkekkxkktx1)()(10010tktktktkeBAecekkkxkkkecty121222212001202)(0222)0(yycBA2022-6-16黄建华制作45例5.2饮酒驾车的数学模型饮酒驾车的数学模型v模型讨论： 当酒是在较短时间内喝时此时有x(0)=x0=D0 ，k0=0 ,y(0)=y0=0 :又酒精浓度为酒精量与体液容积之比，即：)()(122101tktkeekkDkty)()()()()(1212021010tktktktkeeAeevkkDkvtytc2022-6-16黄建华

43、制作46例5.2饮酒驾车的数学模型饮酒驾车的数学模型显然K1K2。当t比较大时, 可认为：利用数表 ：tkAetc2)(tKAtc2ln)(ln时间(小时)0.250.50.7511.522.533.544.55酒精含量153437.5414138.534342925.52520.5时间(小时)678910111213141516酒精含量1917.51412.597.5653.53.522022-6-16黄建华制作47例5.2饮酒驾车的数学模型饮酒驾车的数学模型通过Matlab进行曲线拟合，编程： t= 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 11 12 13 1

44、4 15 16;y1=82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4;y2=log(y1);polyfit(t,y2,1);ans = -0.1940 4.7753可得：A=118.5459,k2=0.1940根据查阅资料可知：一瓶啤酒的酒精量一般为640ml，密度为810mg/ml酒精浓度为84.5%所以两瓶啤酒的酒精总量 由于体重为70kg,体重的65%左右，体液密度为1.05mg/ml，所以可得体液的总体积为2022-6-16黄建华制作48例5.2饮酒驾车的数学模型饮酒驾车的数学模型可得短时间内喝下两瓶啤酒时血液中的酒精含量与时

45、间的关系式如下:)/(33.43310005.110%657030百毫升毫克v114.2)(102101kvkkDkA可得：由5459.118)(114. 21940. 0tteetc2022-6-16黄建华制作49例5.2饮酒驾车的数学模型饮酒驾车的数学模型用Matlab软件画出图形，编程为： t=0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 11 12 13 14 15 16;y1=30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4;a=118.5459;

46、k1=2.114;k2=0.1940;x=0:20;y=a*(exp(-k2*t)-exp(-k1*t);plot(t,y1,+,t,y,r) 2022-6-16黄建华制作50例5.2饮酒驾车的数学模型饮酒驾车的数学模型图形为：2022-6-16黄建华制作51例5.2饮酒驾车的数学模型饮酒驾车的数学模型问题一解答： 假设大李第一次喝酒是在短时间内喝的，根据所建立模型，所以可求得： 当 t=6时，可以求得 c(t)=18.2778毫克/百毫升，小于国家规定的新标准，所以第一次遭遇检查时没有被认定为是饮酒驾驶，见图：27295.59)(114. 21940. 0tteetc百毫升/2778.18)

47、(mgtc2022-6-16黄建华制作52例5.2饮酒驾车的数学模型饮酒驾车的数学模型图形为：2022-6-16黄建华制作53例5.2饮酒驾车的数学模型饮酒驾车的数学模型大李在吃晚饭时又喝了一瓶啤酒，喝酒后体内酒精叠加.当t=14时， c(t)=20.3618毫克/百毫升大于国家新规定的20mg/100ml,所以酒精超标 ：百毫升/2778.18)(mgtc2022-6-16黄建华制作54例5.2饮酒驾车的数学模型饮酒驾车的数学模型问题三解答：短时间内喝酒时，根据所建立模型可知：当 c(t)的导数等于0时，可解得：因为T只与k1,k2有关，从表达式可知当在较短时间内喝酒时浓度达到最大值的时间与喝酒量无关。)()(1202101tktkeeVKKDktc23. 1212ln1lnkkkkT2022-6-16黄建华制作55例5.2饮酒驾车的数学模型饮酒驾车的数学模型其余解答请参考其他资料