第4章程序设计在前面我们已经看到,MATLAB不但可以在命令窗直接输入命令并运行,而且还可以生成自己的程序文件,这就是我们通常说的一类以M为后缀的M文件,本章我们就来研究这类文件的形成方法。
M文件可分分为两大类,一是命令式M文件(也称为脚本文件,script),二是函数式M 文件(function)。
两类文件的区别在于:(1)命令式文件可以直接运行,函数式文件不能直接运行,只能调用。
(2)命令式文件运行时没有输入输出参量,函数式文件在调用时需要进行输入输出参量设置。
(3)命令式文件运行中可以调用工作空间的数据,运行中产生的所有变量为全局变量。
(4)函数式文件不能调用工作空间的数据,运行中产生的所有变量为局部变量。
命令式文件运行中产生的所有变量为全局变量,可以调用和存储到工作空间的数据。
4.1 MATLAB的程序文件-M文件4.1.1 脚本文件(Scripts)当我们需要在命令窗进行大量的命令集合运行时,直接从命令窗口输入比较麻烦,这时就可以将这些命令集合存放在一个脚本文件(Scripts)中,运行时只需要输入其文件名就可以自动执行这些命令集合。
需要注意的是,脚本文件运行所产生的变量都驻留在MATLAB 的工作空间中,同时脚本文件也可以调用工作空间中的数据。
因此,脚本文件所涉及的变量是全局变量。
前几章所涉及到的M文件都是这类脚本文件。
编辑一个脚本文件可以直接在命令窗口的左上角打开编辑窗进行编辑。
4.1.2 函数文件(function)函数式文件(function)的构成(1)函数定义行:Function [输出参量]=gauss(输入参量)(2):完成函数的功能。
(3)函数说明。
(4)函数行注。
从上面构成的情况看,函数式文件实际上是完成输入参量与输出参量的转换,这样的转换是由函数文件名为gauss的文件来完成的。
函数体的功能必须说明清楚输入参量与输出参量的关系。
函数说明是用来解释该函数的功能的,函数行注是对程序行进行说明的。
上面(1)和(2)是必须的。
【例4-1】分析下面函数文件。
%一个数列,任意项等于前两项之和,输入项数可以给出这个数列function [a]=sul(n)if n==1a=1;else if n==2a=2;endb(1)=1;b(2)=2;for i=3:nb(i)=b(i-2)+b(i-1);enda=b;end该函数文件的文件名为sul.m,在第一行给出了该函数的功能,即输入项数就可以自动给出一个满足条件的数列。
定义行规定了输入参数是该数列的项数,输出参数是该数列。
从第三行起,是该函数的主体,主要说明了输入参数与输出参数的关系。
当我们在命令窗或脚本文件中调用该函数就会有结果,如>> p=sul(5)p =1 2 3 5 8【例4-2】分析下面函数文件。
%一个数列的通项公式为a(n+1)=a(n)+n^2,给定任意项的值,求这个数列的后10项,并画离散函数图function [b]=shulie(n,zhi)%zhi为初值b(1)=zhifor m=1:n-1b(1+m)=b(m)+m^2endfor j=n:n+9b(1+j)=b(j)+j^2endx=n:(n+9)stem(x,b(n:n+9))end如果在命令窗输入shulie(5,10)结果为b =10x =5 6 7 8 9 10 11 12 13 14 ans =Columns 1 through 510 11 15 24 40 Columns 6 through 1065 101 150 214 295 Columns 11 through 15395 516 660 829 1025以及4.2 程序的流程控制4.2.1 循环控制1.硬循环语句(for-end)所谓硬循环是指无条件的循环,其结构为for(循环变量)。
end【例4-3】用硬循环语句形成一个5阶对角线元素为2的矩阵。
clearn=5; %每个for和if都必须对应一个endfor i=1:n % i是循环变量for j=1:nif i==ja(i,j)=2;elsea(i,j)=0;endendenda结果为a =2 0 0 0 00 2 0 0 00 0 2 0 00 0 0 2 00 0 0 0 2根据结果,如果要两条主对角线元素都为2,如何调此程序?注意,硬循环语句也可用break来跳出循环。
4.2.2 条件控制1.条件分支语句(if-else-end)如果-那么-否则结构If 条件1语句1elseif 条件2语句2..else语句end注意这一结构的条件优先问题,当有多个条件时,如果条件1不满足,再判断elseif 后面的条件2…,如果所有的条件都不满足,在则执行else 后面的语句。
【例4-5】分析下面条件语句。
x=input('x=?')if ((x+3)<2)y=x*2;elseif (x<2)y=x^2;else(x<1)y=sqrt(x);endy运行结果x=?1.5x =1.5000y =2.25002.条件循环语句(while-end)当条件满足时执行下一语句,否则就到end返回while。
其结构为while(表达式)。
end【例4-4】求阶乘大于或等于99^99的最小整数。
clearn=1;while prod(1:n)<99^99; % prod()向量元素的乘积n=n+1;endn结果为n =1204.2.3 开关与试探控制分支开关语句(switch-case-otherwise-end)结构 switch (开关量),相当于满足什么条件做什么事。
【例4-5】根据开关量的情况来决定计算机计算机执行那个语句。
a=input('a=?')switch acase 1disp('It is raning')case 0disp('It do not konw')case -1disp('It is not raning')otherwisedisp('It is raning?')end4.3 程序设计4.3.1 程序的调试1. 程序的错误一般分为语法和逻辑两种错型。
语法错误时程序无法运行,并在命令窗给予提示。
而逻辑错误只能根据计算的具体情况来判断。
2.根据出错信息调试根据命令窗的提示,注意一般情况不加“;”号调试。
3.设置断点来判断dbstop 格式:dbstop in M文件名 at 行号(也可用鼠标)4.利用keyboard命令在文件中来判断当出现k>> 时return。
5.变量的鼠标观测法当程序代码运行后,每个变量的数据都可以用鼠标放在相应的位置处去观察该变量的数据情况,以判断程序运行的具体执行情况。
6.代码运行的计时方法(1)整段程序代码的计时tic …toc 表示计算tic与toc之间的时间【例4-6】计算机程序运行时间计算。
tica=rand(300);inv(a);%计算逆矩阵toc%注:程序中不需要显示结果的就不显示,这可以节约机时。
结果为>>elapsed_time =0.3900(2)也可以用 etime(t1,t2)来计算t1,t2之间的时间差来完成上述功能。
【例4-7】计算机程序运行时间计算。
t0=clocka=rand(300);inv(a);%计算逆矩阵elapsed_time=etime(clock,t0)%注:程序中不需要显示结果的就不显示,这可以节约机时。
结果仍为>>elapsed_time =0.3900(3)也可以用cputime变量来完成上述功能。
【例4-8】计算机程序运行时间计算。
t0=cputimea=rand(300);inv(a);%计算逆矩阵elapsed_time=cputime-t0%注:程序中不需要显示结果的就不显示,这可以节约机时。
结果为t0 =9.5787e+003elapsed_time =0.27104.3.2 M文件的性能优化1.程序代码的向量化【例4-8】比较下面两个程序段,观察其运行的时间。
程序1t0=cputimen=100000;total=0;for i=1:ntotal=total+1/iendtotalt1=cputime-t0运行结果为total =12.0901t1 =3.1950程序2ticn=100000;a=1:n;total=sum(1./a)toc运行结果为total =12.0901elapsed_time =0.0200发现,程序1和程序2运行结果完全相同,但运行时间程序1是程序2的160倍!2.用矩阵结构进行运算一般情况下,完全采用矩阵运行的方式,MATLAB 的程序与C语言的运行基本相同。
这必须对矩阵非常熟练,例如x=[1 2 3;1 2 1]a=[4 5 6]如果希望将a中的每一个元素乘以x的每一列,用diag(x)。
3.阵的预先配置在一些大的程序运行时,为了变量不至于让计算机不断向内存要空间位置而耽误时间,运行前需要给变量留出一定的空间。
【例4-9】比较下面两个程序段,观察其运行的时间。
程序1tica=[1 2 3;4 5 6;7 8 9];for i=1:100y(i)=det(a^i)endtoc结果为elapsed_time =0.1100程序2tica=[1 2 3;4 5 6;7 8 9];y=zeros(1,100);%分配内存for i=1:100y(i)=det(a^i); endtoc结果为elapsed_time =0.0100练习 41.已知函数13+-=ax x y请研究参数a 在[-1:0.1:1]内函数的根(实部与虚部)随a 的变化关系。
(用图形表示)2.编制一个函数,函数的输入参数为一个任意矩阵或向量,输出参数为该矩阵中不相同的元素个数。
3.编制一个函数文件,函数的输入参数为一个任意矢性场y x e y x g e y x f ρρ),(),(+=ϕ和这个矢性场自变量的取值范围,要求这个函数能自动画出矢性场的平面场分布。
4.设计一个函数文件,要求输入一个任意方阵,该函数能统计出方阵中大于零、等于零和小于零的元素个数。
5.设计一个能实现下面功能的函数文件:输入两个参数:(1)任意2╳n 矩阵(n>8,且第一行为自变量,第二行为对应的函数值);(2)介于任意两自变量之间的值。
输出为该自变量对应的函数值。
(用9次多项式拟和)6. 一个MATLAB 的递归函数fibo.m 来计算fibo 数列,定义如下:fibo (n+2)=(fibo (n+1)+fibo (n ))/2此数列的初条件为fibo (1)=0, fibo (2 )=1n 的最大数为100,要求:(1) 保存你的fibo.m 文件,当在命令窗调用fibo 函数时,不论输入任何整数有正确的输出。