MATLAB数学建模与仿真（第2版·微课视频版）程序代码 第5张 M文件编程

【例5-1】简单函数M文件实例。本例以一个求n的阶乘的函数M文件为例，简单介绍M文件的基本单元。代码如下：functionf=fact(n)%函数定义行，脚本式M文件无此行%Computeafactiorrialvalue.%H1行%FACT(N)returnsthefactionalofN,%Help文本%usuallydenotedbyN!%Putsimple,Fact(N)isPROD(1:N).%注释f=prod(1:n);%函数体或脚本主体【例5-2】Help文本查看。本例演示通过help命令查看M文件中的帮助文本的过程。>>helpfactComputeafactiorrialvalue.%H1行fact(N)returnsthefactionalofN,%Help文本usuallydenotedbyN!Putsimple,Fact(N)isPROD(1:N).%注释以上命令显示fact函数文件的注释行，直到第1个非注释行（即空行）结束。键入lookfor命令可见：>>lookforfactfact-Computeafactiorrialvalue.ipjfact-Hankelmatrixwithfactorialelements.chol-Choleskyfactorization.cholupdate-Rank1updatetoCholeskyfactorization.ldl-BlockLDL'factorizationforHermitianindefinitematrices.lu-LUfactorization…【例5-3】注释块操作符示例。将例5-1的fact函数中的多行注释改写为注释块。functionf=fact(n)%函数定义行，脚本式M文件无此行 %{Computeafactiorrialvalue.%H1行 FACT(N)returnsthefactionalofN,%Help文本 usuallydenotedbyN! Putsimple,Fact(N)isPROD(1:N).%注释%} f=prod(1:n); %函数体或脚本主体【例5-4】编写求平均值与标准差的脚本文件stat1.m和函数文件stat2.m。单击MATLAB命令窗口工具栏中的按钮，打开M文件编辑器。在两个文件各自的M文件编辑器中分别输入以下程序代码：%stat1.m脚本文件%求矩阵x每列的平均值和标准差[m,n]=size(x);ifm==1m=n;ends1=sum(x);s2=sum(x.^2);mean1=s1/m;stdev=sqrt(s2/m-mean1.^2);%stat2.m函数文件function[mean1,stdev]=stat2(x)%求矩阵x每列的平均值和标准差[m,n]=size(x);ifm==1m=n;ends1=sum(x);s2=sum(x.^2);mean1=s1/m;stdev=sqrt(s2/m-mean1.^2);在命令窗口依次输入如下命令：>>clear>>x=rand(4,4)+2;>>stat1%执行stat1.m后，观察基本空间中的变量情况。>>whos%可见脚本文件所产生的所有变量都返回工作空间NameSizeBytesClassAttributesm1x18doublemean11x432doublen1x18doubles11x432doubles21x432doublestdev1x432doublex4x4128double>>disp([mean1;stdev])%观察计算结果2.43512.53232.57282.59140.33160.25770.10850.1857在命令窗口依次输入如下命令：>>clearmns1s2mean1stdev>>[m1,std1]=stat2(x);%执行stat2.m后，观察基本基本空间中的变量情况>>whos%只增加了由函数返回的结果NameSizeBytesClassm11x432doublestd11x432doublex4x4128double>>disp([m1;std1])%观察计算结果，和stat1.m一致2.43512.53232.57282.59140.33160.25770.10850.1857【例5-5】将文件MATLAB.mat中的变量导入到工作区中。（1）新建一个MATLAB.mat文件。打开MATALB，在选项卡中单击图标，即可打开一个M文件编辑器，输入以下程序代码。保存为文件Matlab.m，并单击。%生成基础测量数据x=-3*pi:3*pi;y=x;[X,Y]=meshgrid(x,y);R=sqrt(X.^2+Y.^2)+eps;Z=sin(R)./R;[dzdx,dzdy]=gradient(Z);dzdr=sqrt(dzdx.^2+dzdy.^2);%绘制基础数据图形surf(X,Y,Z,abs(dzdr))colormap(spring)alphamap('rampup')colorbar%进行二维插值运算xi=linspace(-3*pi,3*pi,100);yi=linspace(-3*pi,3*pi,100);[XI,YI]=meshgrid(xi,yi);ZI=interp2(X,Y,Z,XI,YI,'cubic');%绘制插值后的数据图形figuresurf(XI,YI,ZI)colormap(spring)alphamap('rampup')colorbarsaveMATLAB（2）用命令whos–file查看该文件中的内容：>>whos-fileMATLAB.matNameSizeBytesClassAttributesR19x192888doubleX19x192888doubleXI100x10080000doubleY19x192888doubleYI100x10080000doubleZ19x192888doubleZI100x10080000doubledzdr19x192888doubledzdx19x192888doubledzdy19x192888doublex1x19152doublexi1x100800doubley1x19152doubleyi1x100800double（3）将该文件中的变量导入到工作区中：>>clear>>loadMATLAB.mat（4）该命令执行后，可以在工作区窗口中看见这些变量，如图5-6所示。图5-6导入变量后的工作区视图（5）接下来用户可以访问这些变量：>>dzdydzdy=列1至11-0.0459-0.0568-0.0418-0.00650.03530.07060.09220.10080.10140.10010.1002…【例5-6】从文件中导入数据。>>imported_data=importdata('MATLAB.mat')import_data=包含以下字段的struct:R:[19×19double]X:[19×19double]…y:[1×19double]yi:[1×100double]【例5-7】在MATLAB中使用fopen命令打开磁盘文件。以读写的方式打开磁盘文件fgetl.m。在MATLAB的命令窗口中输入以下代码：>>[fid,message]=fopen('fgetl.m','r+')查看程序结果。在输入以上代码后，得到的结果如下：fid=-1message= Nosuchfileordirectory【例5-8】open与load的比较。>>clear>>A=magic(3);>>B=rand(3);>>save正在保存到:D:\Document\MATLAB\MATLAB.mat>>clear>>load('MATLAB.mat')>>AA= 816 357 492>>BB=0.25750.81430.35000.84070.24350.19660.25430.92930.2511>>clear>>open('MATLAB.mat')ans=包含以下字段的struct:A:[3×3double]B:[3×3double]>>struc1=ans;>>struct1=ans>>struct1.Aans= 816 357492>>struct1.Bans=0.25750.81430.35000.84070.24350.19660.25430.92930.2511【例5-9】在MATLAB中关闭对应的磁盘文件。创建文件fgetl.m，然后删除该文件。在MATLAB的命令窗口中输入以下代码：>>[fid,message]=fopen('fget1.m','w')>>deletefgetl.m查看程序代码的结果。以上程序代码可以得到如下结果：警告：未找到文件或者权限被拒绝以上结果表明，当用户使用fopen命令创建了对应的空白文件fgetl.m，并打开对应的文件后，如果在关闭该文件前试图删除文件，系统会提示用户删除命令被拒绝。首先关闭文件。然后删除该文件。在MATLAB的命令窗口中输人以下代码：>>status=fclose(fid);>>deletefgetl.m;>>[fid,message]=fopen('fget1.m','r+');查看程序代码的结果。当用户输入以上程序代码后，得到的结果如下：fid=-1message= Nosuchfileordirectory【例5-10】显示basic_matrix.txt的内容，然后使用readmatrix将数据导入矩阵。在命令窗口中输入：>>typebasic_matrix.txt6,8,3,15,4,7,31,6,7,104,2,8,22,7,5,9>>M=readmatrix('basic_matrix.txt')M=5×4683154731671042822759 从以上结果可以看出，使用readmatrix命令可以读取文件中的数据，并且将读取的数据导入矩阵。【例5-11】创建一个矩阵，将其写入以逗号分隔的文本文件，然后用不同分隔符将该矩阵写入另一个文本文件。在工作区中创建一个矩阵。>>M=magic(5)M=17241815235714164613202210121921311182529 将矩阵写入逗号分隔的文本文件，并显示文件内容。writematrix函数将输出名为M.txt的文本文件。>>writematrix(M)>>type'M.txt'17,24,1,8,1523,5,7,14,164,6,13,20,2210,12,19,21,311,18,25,2,9要用不同分隔符将同一矩阵写入文本文件，请使用'Delimiter'名称-值对组。>>writematrix(M,'M_tab.txt','Delimiter','tab')>>type'M_tab.txt'1724 1 8 1523 5 7 14 164 6 13 20 2210 1219 21 311 1825 2 9【例5-12】在MATLAB中使用textscan命令读取文本文件。查看原始的数据文件。在本实例中，用户需要读取的文件是textscan.dat，其文件中包含的数据如下：SallyLevel112.34451.23e10infNaNYesJoeLevel223.54609e19-inf0.001NoBillLevel334.90122e510100No使用命令读取该数据文件。在MATLAB的命令窗口中输入以下程序代码：>>fid=fopen('textscan.dat');>>Cl=textscan(fid,'%s%s%f32%d8%u%f%f%s');>>fclose(fid);查看程序代码的结果。在命令窗口中输入变量名称，得到的结果如下：>>whosClNameSizeBytesClassAttributesCl1x81917cell可以看出，使用textscan命令得到的结果将会被存储在元胞数组中，该元胞数组包含的列数就是原始数据文件中使用分隔符隔开的数据列。查看C1数组的结果。在命令窗口中输入以下程序代码：>>fori=1:8disp(C1{i}');end查看程序代码的结果。输入代码后，按Enter键，得到的结果如下：{'Sally'}{'Joe'}{'Bill'}{'Level1'}{'Level2'}{'Level3'}12.340023.540034.900045601242949672954294967295200000Inf-Inf10NaN0.0010100.0000{'Yes'}{'No'}{'No'}读取原始文件，并忽略第3列数据。在命令窗口中输入以下程序代码：>>fid=fopen('textscan.dat');>>C2=textscan(fid,'%7c%6s%*f%d8%u%f%f%s');>>fclose(fid);查看C2的属性。在命令窗口中输入以下程序代码：>>whosC2NameSizeBytesClassAttributesC21x71801cell查看C2数组的结果。在命令窗口中输入以下程序代码：>>fori=1:7disp(C2{i});end查看程序代码的结果。输入代码后，按Enter键，得到的结果如下：{'Sally'}{'Joe'}{'Bill'}{'Level1'}{'Level2'}{'Level3'}45601242949672954294967295200000Inf-Inf10NaN0.0010100.0000{'Yes'}{'No'}{'No'}'从以上结果可以看出，当在textscan命令中使用textscan(fid,'%7c%6s%*f%d8%u%f%f%s')后，其中“%*f”所替代的对应数据列会被跳过，不被读入。仅读取原始文件的第一列数据。在命令窗口中输入以下程序代码：>>fid=fopen('textscan.dat');>>names=textscan(fid,'%s%*[^\n]');>>fclose(fid);查看names的属性。在命令窗口中输入以下程序代码：>>whosnamesNameSizeBytesClassAttributesnames1x1440cell查看程序代码的结果。在命令窗口中输入以下程序代码：>>B=names{1}B=3×1cell数组{'Sally'}{'Joe'}{'Bill'}消除原始第3列数据前面的标签。在命令窗口中输入以下程序代码：>>fid=fopen('textscan.dat');>>C3=textscan(fid,'%sLevel%u8%f32%d8%u%f%f%s');>>fclose(fid);查看程序代码的结果。输入代码后，按Enter键，得到的结果如下：>>whosC3NameSizeBytesClassAttributesC31x81572cell查看C3的属性。在命令窗口中输入以下程序代码：>>fori=1:8disp(C3{i});end查看程序代码的结果。输入代码后，按Enter键，得到的结果如下：{'Sally'}{'Joe'}{'Bill'}12312.340023.540034.900045601242949672954294967295200000Inf-Inf10NaN0.0010100.0000{'Yes'}{'No'}{'No'}'从以上结果中可以看出，相对于原始的第二列数据，通过该命令得到的第二列数据清除了字符串Level，只留下了数值代码。【例5-13】顺序结构示例。实现计算a与b的和与积，相乘后减去c的功能。编写M文件Ex_5_16.m，代码如下所示。a=1;b=2;c=3;ans1=a+bans2=a*bans3=ans2*ans1-c单击F5快捷键或“运行”按钮运行该文件，或者将其在当前目录下保存为Ex_5_16.m，然后在命窗口中输入Ex_5_16.m并运行，得到如下结果：>>Ex_5_16ans1=3ans2=2ans3=3【例5-14】顺序结构示例。在MATLAB中，使用顺序结构编写绘制函数的图形。将文件保存为Ex_5_17.m。%定义符号变量t和taosymsttao%定义积分表达式y=exp(-t/3)*cos(1/2*t);%对表达式进行积分S=subs(int(y,0,tao),tao,t);%绘制积分图形ezplot(s,[0,4*pi]);grid返回到MATLAB的命令窗口，输入Ex_5_17，然后按Enter键，得到的结果如图5-7所示。图5-7得到的程序结果在上面的程序代码中，首先定义符号变量，然后定义积分表达式，进行积分运算，最后调用ezplot命令绘制积分函数的图形。这样的程序代码流程符合逻辑顺序，而且容易阅读，容易理解，这是顺序结构的重要优点。【例5-15】if…end语句的实例。文件保存为EX_5_18.m。a=6;ifrem(a,2)==0%判断a是否是偶数disp('aiseven')b=a/2end本例中的程序首先判断a是否是偶数，因为a的值为6，所以命令rem(a,2)==0返回逻辑值true。然后程序运行if语句之内的程序段,得出如下结果：>>Ex_5_18aisevenb=3【例5-16】if…else…end语句使用示例。 文件保存为Ex_5_19.m。ifa>bdisp('aisbiggerthenb')%a>b则执行此句y=a;elsedisp('aisnotbiggerthenb')%a<=b则执行此句y=b;end【例5-17】if…elseif…else…end语句使用示例。ifn<0%如果n是负数，则显示错误信息disp('Inputmustbepositive');elseifrem(n,2)==0%如果n是偶数，则除以2A=n/2;elseA=(n+1)/2;%如果n是奇数，则加1，然后除以2end【例5-18】在MATLAB中，使用if分支结构编写求解一元二次方程ax首先分析分支结构的判断条件。根据基础数学知识可知，一元二次方程ax2+bx+c=0的根的性质直接取决于判别式∆=b2-4ac的数值。当∆单击MATLAB命令窗口工具栏中的按钮，打开M文件编辑器。在M文件编辑器中输入以下程序代码：%scriptfilecale_root.m%purpose:%Thisprogramsolvesfortherootsofaquadraticequation%oftheforma*x^2+b*x+c=0.Itcalcaulatestheanswersof%rootstheequationposseses.%Definevariables:%a=coefficientofx^2%bcoefficientofx%cconstantterm%x1firstrootoftheequation%x2secondrootoftheequationdisp('Thisprogramsolvesfortherootsofaquadraticequation');disp('oftheforma*x^2+b*x+c=0');a=input('EnterthecoefficientA:');b=input('EnterthecoefficientB:');c=input('EnterthecoefficientC:');discriminant=b^2-4*a*c;%如果判别式大于0%则根据二元方程的公式得出两个不同的实数解ifdiscriminant>0x1=(-b+sqrt(discriminant))/(2*a);x2=(-b-sqrt(discriminant))/(2*a);%在命令窗口显示求解结果disp('Thisequationhastworealroots');fprintf('x1=%f\n',x1);fprintf('x2=%f\n',x2);%当判别式等于0，则返回两个相同的实数根elseifdiscriminant==0x1=-b/(2*a);disp('Thisequationhastwoidenticalroots');fprintf('x1=x2=%f\n',x1);%当判别式小于0，则返回两个须根elsereal_part=-b/(2*a);image_part=sqrt(abs(discriminant))/(2*a);disp('Thisequationhastwocomplexroots');fprintf('x1=%f+i%f\n',real_part,image_part);fprintf('x2=%f-i%f\n',real_part,image_part);end单击M文件编辑器中的“保存”按钮，将以上程序代码保存为calc_root.m。返回到MATLAB的命令窗口，输入calc_root，然后按Enter键，根据程序代码的提示，依次输入方程的系数，得到的结果如下：>>calc_rootThisprogramsolvesfortherootsofaquadraticequationoftheforma*x^2+b*x+c=0EnterthecoefficientA:1EnterthecoefficientB:5EnterthecoefficientC:6Thisequationhastworealrootsx1=-2.000000x2=-3.000000>>calc_rootThisprogramsolvesfortherootsofaquadraticequationoftheforma*x^2+b*x+c=0EnterthecoefficientA:2EnterthecoefficientB:3EnterthecoefficientC:4Thisequationhastwocomplexrootsx1=-0.750000+i1.198958x2=-0.750000-i1.198958>>calc_rootThisprogramsolvesfortherootsofaquadraticequationoftheforma*x^2+b*x+c=0EnterthecoefficientA:1EnterthecoefficientB:4EnterthecoefficientC:4Thisequationhastwoidenticalrootsx1=x2=-2.000000【例5-19】switch…case…end示例1。switchvarcase1%判断var是不是1disp('1')case{2,3,4}%判断var是不是2,3,4disp('2or3or4')case5%判断var是不是5disp('5')otherwise%其他情况disp('somethingelse')end【例5-20】switch…case…end示例2。文件保存为Ex_5_20.m。clear;%划分区域：满分(100)，优秀(90-99)，良好(80-89)，及格(60-79)，不及格(<60)fori=1:10;a{i}=89+i;b{i}=79+i;c{i}=69+i;d{i}=59+i;end;c=[d,c];Name={'Jack','Marry','Peter','Rose','Tom'};Mark={72,83,56,94,100};Rank=cell(1,5);%3个数组，且都是(1x5)维的%创建一个含5个元素的构架数组S，它有三个域S=struct('Name',Name,'Marks',Mark,'Rank',Rank);%根据学生的分数，求出相应的等级fori=1:5switchS(i).Markscase100S(i).Rank='满分';caseaS(i).Rank='优秀';casebS(i).Rank='良好';casecS(i).Rank='及格';otherwiseS(i).Rank='不及格';endend%将学生姓名，得分，等级等信息打印出来disp(['学生姓名','得分','等级']);disp('')fori=1:5;disp([S(i).Name,blanks(6),num2str(S(i).Marks),blanks(6),S(i).Rank]);end;运行结果为：>>Ex_5_23学生姓名得分等级Jack72及格Marry83良好Peter56不及格Rose94优秀Tom100满分【例5-21】try…catch结构应用实例。对3×3魔方阵的行进行引用，当“行下标”超出魔方阵的最大行数时，将改向对最后一行的援引，并显示“出错”警告。 单击MATLAB命令窗口工具栏中的按钮，打开M文件编辑器。在M文件编辑器中输入以下程序代码：clearN=4;A=magic(3);%设置3行3列矩阵Atry A_N=A(N,:)%取A的第N行元素catchA_end=A(end,:)%如果取A(N,:)出错，则改取A的最后一行endlasterr%显示出错原因单击M文件编辑器中的“保存”按钮，将以上程序代码保存为“Ex_5_21.m”。返回到MATLAB的命令窗口，输入Ex_5_21，然后按Enter键，根据程序代码的提示，依次输入方程的系数，得到的结果如下：>>Ex_5_21A_end=492ans='位置1处的索引超出数组边界(不能超出3)。'【例5-22】for循环使用示例1。文件保存为EX_5_22.m。x=ones(1,6)forn=2:6%循环控制x(n)=2*x(n–1)%循环体end运行后可得到如下结果：>>Ex_5_22x=111111x=121111x=124111x=124811x=1248161x=12481632【例5-23】for循环使用示例2。文件保存为EX_5_23.m。form=1:5forn=1:10A(m,n)=1/(m+n-1);%使用循环体给变量A赋值endend运行后可得到如下结果：>>Ex_5_23>>AA=1.00000.50000.33330.25000.20000.16670.14290.12500.11110.10000.50000.33330.25000.20000.16670.14290.12500.11110.10000.09090.33330.25000.20000.16670.14290.12500.11110.10000.09090.08330.25000.20000.16670.14290.12500.11110.10000.09090.08330.07690.20000.16670.14290.12500.11110.10000.09090.08330.07690.0714【例5-24】while循环使用示例。文件保存为Ex_5_24.m。i=1;whilei<10%i小于10进行循环x(i)=i^3;%循环体内的计算i=i+1;%表达式值的改变end运行后可以得到如下结果:>>Ex_5_24>>xx=182764125216343512729>>ii=10【例5-25】多种循环体的嵌套使用示例。文件保存为Ex_5_25.m。clearclcfori=1:1:6%行号循环，从1到6j=6;whilej>0%列号循环，从6到1x(i,j)=i-j;%矩阵x的第i行第j列元素值为其行号的差ifx(i,j)<0x(i,j)=-x(i,j);%当x(i,j)为负数时，取其相反数endj=j-1;endend运行Ex_5_25.m文件，可以得到如下结果：>>Ex_5_25>>xx=012345101234210123321012432101543210【例5-26】continue命令的使用方法。单击MATLAB命令窗口工具栏中的按钮，打开M文件编辑器。在M文件编辑器中输入以下程序代码：forii=1:9%ifii==3之所以在这段程序代码中的第二段代码前面添%continue加了%，是为了首先将其当作注释，不运行这段代码，%end在后面的程序中将注释符号删除后就可以重新运行。fprintf('ii=%d\n',ii);ifii==5breakendenddisp('Theendofloop');将以上代码保存为break_continue.m文件，然后在MATLAB的命令窗口中输入break_continue，按Enter键，就可以得到对应的结果：>>break_continueii=1ii=2ii=3ii=4ii=5Theendofloop打开break_continue.m文件，然后在编辑器中修改其代码，得到的结果如下：forii=1:9ifii==3 continueendfprintf('ii=%d\n',ii);ifii==5breakendenddisp('Theendofloop');在MATLAB的命令窗口中输入“break_continue”，按“Enter”键，就可以得到对应的结果：>>break_continueii=1ii=2ii=4ii=5Theendofloop注意：在上面程序代码中使用了break语句，其功能就是跳出相应的程序代码。【例5-27】在MATLAB中寻求Fibonacci数组中第一个大于700的元素以及其数组标号。单击MATLAB命令窗口工具栏中的按钮，打开M文件编辑器。在M文件编辑器中输入以下程序代码：n=50;a=ones(1,n);fori=3:na(i)=a(i–1)+a(i–2);ifa(i)>=700a(i)break;endendi将以上代码保存为Fib.m文件，然后在MATLAB的命令窗口中输入Fib，按Enter键，就可以得到对应的结果：>>Fibans=987i=16从以上结果可以看出，在Fibonacci数组中第一个大于700的数值是987，其对应的数组标号是16。【例5-28】在MATLAB中演示如何使用input函数。 单击MATLAB命令窗口工具栏中的按钮，打开M文件编辑器。在M文件编辑器中输入以下程序代码：functiontest_input()%在以上程序代码中，使用isempty来接收用户输入的Enter键%当什么字符都不输入的时候，默认认为用户输入的是Y。reply=input('Doyouwantmore?Y/N[Y]:','s');ifisempty(reply)reply='Y';endifreply=='Y'disp('youhaveselectedmoreinformation');elsedisp('youhaveselectedtheend');end将以上代码保存为test_input.m文件。在MATLAB的命令窗口中输入test_input，然后按Enter键，就可以得到对应的结果：>>test_inputDoyouwantmore?Y/N[Y]:youhaveselectedmoreinformation>>test_inputDoyouwantmore?Y/N[Y]:Yyouhaveselectedmoreinformation>>test_inputDoyouwantmore?Y/N[Y]:Nyouhaveselectedtheend【例5-29】在MATLAB中，演示如何使用keyboard命令。 在MATLAB的命令窗口中输入以下内容：>>keyboardK>>forii=1:9ifii==3continueendfprintf('ii=%d\n',ii);ifii==5breakendendii=1ii=2ii=4ii=5K>>dbcont>> 从以上程序代码可以看出，当输入keyboard命令后，在提示符的前面会显示k提示符，而当输入dbcont后，提示符恢复正常。【例5-30】使用不同的警告样式，查看MATLAB的不同错误提示模式。单击MATLAB命令窗口工具栏中的按钮，打开M文件编辑器。在M文件编辑器中输入以下程序代码：文件保存为error_message.m。%Scriptfileerror_message.m%%Purpose:%Tocalculatemeanandthestandarddeviationof%aninputdatasetcontainingandarbitrarynumber%ofinputvalues.%%Definevariables:%nThenumberofinputsamples%std_devThestandarddevationoftheinputsamples%sum1Thesumoftheinputvalues%sum2Thesumofthesquaresoftheinputvalues%xinputdatavalue%xvarTheaverageoftheinputsamples%Initalizevariablessum1=0;sum2=0;%Getthenumberofpointstoinputn=input('Enterthenumberofpoints:'); %Checktoseeifwehaveenoughinputdata.ifn<2errordlg('Notenoughinputdata');else%constructforloopforii=1:nx=input('Entervalue:');sum1=sum1+x;sum2=sum2+x^2;end %Calculatethemeanandstandarddeviation xvar=sum1/n; std_dev=sqrt((n*sum2-sum1^2)/(n*(n-1))); %Printtheresult fprintf('Themeanofthisdatasetis:%f\n',xvar); fprintf('Thestandarddeviationis:%f\n',std_dev); fprintf('Thenumberofdatais:%d\n',n);end返回到MATLAB的命令窗口，输入error_message，然后输入数值1，得到的结果如图5-8所示。当输入的数值总数小于2时，MATLAB调用错误信息对话框。单击对话框中的OK按钮后，将会自动退出程序代码。图5-8显示错误信息打开error_message.m文件，在编辑器中修改其程序代码，然后保存相应的程序代码，修改的程序代码如下：sum1=0;sum2=0;%Getthenumberofpointstoinputn=input('Enterthenumberofpoints:');%Checktoseeifwehaveenoughinputdata.ifn<2error('Notenoughinputdata');else…end返回到MATLAB的命令窗口中，输入error_message，然后输入数值1，得到如下的结果：>>error_messageEnterthenumberofpoints:1???Errorusing==>error_messageat24Notenoughinputdata打开error_message.m文件，在编辑器中修改其程序代码，然后保存相应的程序代码。修改的程序代码如下：sum1=0;sum2=0;%Getthenumberofpointstoinputn=input('Enterthenumberofpoints:');%Checktoseeifwehaveenoughinputdata.ifn<2warning('Notenoughinputdata');else…end返回到MATLAB的命令窗口，输入error_message，然后输入数值1，得到如下的结果：>>error_messageEnterthenumberofpoints:1Warning:Notenoughinputdata>>Inerror_messageat24【例5-31】编写一个M函数文件，它具有以下功能：①根据指定的半径，画出蓝色圆周线；②可以通过输入字符串，改变圆周线的颜色、线型；③假若需要输出圆面积，则绘出圆。本例演示：M函数文件的典型结构；指令nargin、nargout的使用和函数输入输出参数个数的变化；switch…case控制结构的应用示例；if…elseif…else的应用示例；error的使用。单击MATLAB命令窗口工具栏中的按钮，打开M文件编辑器。在M文件编辑器中输入以下程序代码：function[S,L]=sexangle(N,R,str)%sexangle.mTheareaandperimeterofaregularpolygon% (正多边形面积和周长)%N Thenumberofsides%R Thecireumradius%strAlinespecificationtodeterminelinetype/color%S Theareaoftheregularpolygon%L Theperimetertheregularpolygon%sexangle 用蓝实线画半径为1的圆%sexangle(N) 用蓝实线画外接半径为1的正N边形%sexangle(N,R) 用蓝实线画外接半径为R的正N边形%sexangle(N,R,str)用str指定的线画外接半径为R的正N边形%S=sexangle(…)给出多边形面积S，并画相应正多边形填色图%[S,L]=sexangle(…)给出多边形面积S和周长L,并画相应正多边形填色图switchnargincase0N=100;R=1;str='-b';case1R=1;str='-b';case2str='-b';case3;%不进行任何变量操作，直接跳出switch_case控制结构otherwiseerror('输入量太多。');end;t=0:2*pi/N:2*pi;x=R*sin(t);y=R*cos(t);ifnargout==0plot(x,y,str);elseifnargout>2error('输出量太多。');elseS=N*R*R*sin(2*pi/N)/2;%多边形面积L=2*N*R*sin(pi/N);%多边形的周长fill(x,y,str)endaxisequalsquareboxonshg单击M文件编辑器中的“保存”按钮，将以上程序代码保存为sexangle.m。返回到MATLAB的命令窗口，根据不同的解题条件输入不同的值，然后按Enter键，得到不同的图。例如：sexangle(6,2,'-g')sexangle()sexangle(6)sexangle(6,2)以上命令的执行结果见图5-9至图5-12.图5-9用绿色线画外接半径为2的正6边形图5-10用蓝实线画半径为1的圆图5-11用蓝实线画外接半径为1的正6边形图5-12用蓝实线画外接半径为2的正6边形【例5-32】在MATLAB中，通过程序来演示小球的抛物线轨迹。分析小球的抛物线轨迹模型。假定用户抛小球的速度，也就是小球的初始速度是v0，小球的抛射初始角度是θ。根据基础的物理知识可知，小球在水平和垂直方向上的速度分量分别为在本实例中，程序代码需要求解的是抛物线轨迹上水平距离的最长距离，根据相关知识，其距离的求解公式如下：在以上公式中，g代表的是重力加速度。在本实例中该参数选择的数值为-9.82。而对应的小球在垂直方向上的最高距离为根据本实例的要求，可以输入抛射小球的初始速度，然后得出相应的计算数据。单击MATLAB命令窗口工具栏中的按钮，打开M文件编辑器，输入以下程序代码：%Scriptfileball.m%%Purpose:%Thisprogramcalculatesthedistancetraveledbyaball%throwataspeciifiedangle"theta"andaspecifiedvelocity%v0fromapointonthesurfaceoftheearch,ignoringairfrictionandtheearth’scurvature.%Itcalculatestheangleyeiledingmaximunrange,andalsoplotsselectedtrajectories.%Definevariable:%convdegreestoradiansconvfactor%gravThegravityaccel%ii,jjLoopindex%indexThemaximumrangeinarray%maxangleTheanglethatgivesthemaximumrange%maxrangeMaximumrange%rangerangheforaspecifiedangle%timeTime%thetaInitalangle%fly_timethetotletrajectorytime%voTheinitialvelocity%vxox-componentoftheinitialvelocity%vyoy-componentoftheinitialvelocity%xx-positionofball%yy-positionofball%定义常数数值conv=pi/180;grav=-9.82;vo=input('Entertheinitialvelocity:');range=zeros(1,91);%计算最大的水平距离forii=1:91theta=ii-1;vxo=vo*cos(theta*conv);vyo=vo*sin(theta*conv);max_time=-2*vyo/grav;range(ii)=vxo*max_time;end%显示计算水平距离的列表fprintf('Rangeversusangletheta"\n');forii=1:5:91theta=ii-1;fprintf('%2d%8.4f\n',theta,range(ii));end%计算最大的角度和水平距离[maxrangeindex]=max(range);maxangle=index-1;fprintf('\nMaxrangeis%8.4fat%2ddegress.\n',maxrange,maxangle);%绘制轨迹图形forii=5:10:80theta=ii;vxo=vo*cos(theta*conv);vyo=vo*sin(theta*conv);max_time=-2*vyo/grav;%计算小球轨迹的x,y坐标数值x=zeros(1,21);y=zeros(1,21);forjj=1:21;time=(jj-1)*max_time/20;x(jj)=vxo*time;y(jj)=vyo*time+0.5*grav*time^2;endplot(x,y,'g');ifii==5holdon;endend%添加图形的标题和坐标轴名称title('\bfTrajectoryofBallvsInitalAngle\theta');xlabel('\bf\itx\rm\bf(meters)');ylabel('\bf\ity\rm\bf(meters)');axis([0max(range)+50-vo^2/2/grav]);gridon;%绘制最大水平的轨迹图形vxo=vo*cos(maxangle*conv);vyo=vo*sin(maxangle*conv);max_time=-2*vyo/grav;%Calculatethe(x,y)x=zeros(1,21);y=zeros(1,21);forjj=1:21time=(jj-1)*max_time/20;x(jj)=vxo*time;y(jj)=vyo*time+0.5*grav*time^2;endplot(x,y,'r','Linewidth',2);holdoff单击M文件编辑器中的“保存”按钮，将以上程序代码保存为ball.m。返回MATLAB的命令窗口，输入ball，然后按Enter键，根据程序代码的提示，依次输入不同的值，得到的结果如下：>>ballEntertheinitialvelocity:20Rangeversusangletheta"00.000057.07321013.93161520.36662026.18282531.20343035.27603538.27674040.11444540.73325040.11445538.27676035.27606531.20347026.18287520.36668013.9316857.0732900.0000Maxrangeis40.7332at45degress.除了以上数值结果之外，MATLAB还会绘制相应的图形结果，如图5-13所示。图5-13初始速度为20时的轨迹修改初始速度数值，将其改为45，得到的结果如下：>>ballEntertheinitialvelocity:45Rangeversusangletheta"00.0000535.80831070.528615103.105920132.550425157.967430178.584735193.775740203.079045206.211850203.079055193.775760178.584765157.967470132.550475103.10598070.52868535.8083900.0000Maxrangeis206.2118at45degress.同时，MATLAB会给出对应的图形结果,如图5-14所示。图5-14初始速度为45时的轨迹【例5-33】通过M脚本文件，画出下列分段函数所表示的图形。a=2;b=2;clf;x=-a:0.2:a;y=-b:0.2:b;fori=1:length(y)forj=1:length(x)ifx(j)+y(i)>1z(i,j)=0.5457*exp(-0.75*y(i)^2-3.75*x(j)^2-1.5*x(j));elseifx(j)+y(j)<=1z(i,j)=0.5457*exp(-0.75*y(i)^2-3.75*x(j)^2+1.5*x(j));elsez(i,j)=0.7575*exp(-y(i)^2-6*x(j)^2);endendendaxis([-a,a,-b,b,min(min(z)),max(max(z))]);colormap(flipud(winter));surf(x,y,z);将以上内容的M文件Ex_5_36.m保存在当前目录下，然后在命令行输入该M文件的文件名，Ex_5_36，或者按下F5快捷键，或者打开文件后单击M文件编辑器的“运行”按钮，即可运行该文件，运行结果如图5-15所示。图5-15分段函数所对应的曲面【例5-35】子函数的帮助文本查看示例。>>helpnewsstats>meanCallculateaverage.>>helpmeanmean-数组的均值此MATLAB函数返回A沿大小不等于1的第一个数组维度的元素的均值。M=mean(A)M=mean(A,'all')M=mean(A,dim)M=mean(A,vecdim)M=mean(___,outtype)M=mean(___,nanflag)另请参阅median,mode,std,sum,varmean的文档名为mean的其他函数【例5-35】子函数的帮助文本查看示例。>>helpnewsstats>