第二讲Matlab编程与作图
第一部分Matlab程序设计初步
Matlab除了指令行操作的直接交互外，作为一种高级应用软件还提供了自己的编程语言。通过编写Matlab程序，可以更加方便地调用Matlab提供的各种功能强大的函数库，使得程序能完成复杂的运算处理大量的数值数据。
1、M文件简介
Matlab提供了丰富的编程语言，使得用户可以将一连串的命令写入文件，然后使用简单的函数来执行这些命令。文件被保存为文本文件，后缀为.m，比如说dblquad.m，因此Matlab的程序通常被称为M 文件。
M文件是一个文本文件，可以使用各种文本编辑器对它进行编辑和修改，比如Windows操作系统自带的记事本，也可以用Matlab 内建的M文件编辑器。
M文件分为两类，一类称为脚本（Scripts），类似于批处理文件，相当于将在Matlab命令窗口中执行的一系列指令放在一个文件中，当在命令窗口调用该文件名时，则按顺序执行其中的命令集。
例：编写求10！的程序。
另一类M文件称为函数（Function），它可以接受输入变量，并将运算结果送至输出变量，类似于数学中的函数y=f(x)。
函数M文件的基本结构：
function f=fact(n) 函数定义行
%Compute a factorial value.
%FACT(N) returns the factorial of N, 帮助文档%usually denoted by N!
%Put simply,FACT(N) is PROD(1:N), 注释
f=prod(1:n); 函数体2、运算符
关系运算符：<, <=, >, >=, = =, ~=
逻辑运算符：与（&），或（|），非（~）
例：编写分段函数
21
() 1 -1<1
321
x x
f x x
x x
?>
?
=≤
?
?+≤-
?
%myfun1.m
function y=myfun1(x)
y=(x.^2).*(x>1)+(x>-1& x<=1)+(3+2*x).*(x<=-1);
注意：1.函数名与变量名的命名法则相同，要求以字母开头，后接字母或下划线；2.函数名与保存的文件名最好一致。
3、控制流
所有的计算机编程语言都提供了控制程序流执行程序的语法，Matlab也不例外。所有的控制流语法都以end 结尾。
⑴for 循环语句
语法：for 循环变量=数组
指令组；
end
解释：对于循环变量依次取数组中的值，循环执行指令组直到循环变量遍历数组。数组最常用的形式是初值：步长：终值。 例：构造Hibert 矩阵 ⑵while 循环 语法：while 条件式 指令组； end
解释：当条件式满足，循环执行指令组直到条件式不满足。
使用while 语句要注意避免出现死循环。 例：利用迭代公式11()2k k k
a x x x +=+
精度。[Sqrt.m] ⑶分支语句
语法：if 条件表达式1 指令组1； [ elseif 条件表达式2 指令组2；] [·······] [ else 指令组k;] end
解释：如果条件表达式1满足，则执行指令组1，且结束该语句；否
则检查条件表达式2，若满足则执行指令组2，且结束该语句；······；若所有的条件都不满足，则执行指令组k，并结束该语句。
例：用条件语句编写分段函数[myfun2.m]
⑷开关语句
语法：switch 分支变量
case 值1
指令组1；
case 值2
指令组2；
··········
otherwise
指令组k;
end
解释:若分支变量的值取值1，则执行指令组1，且结束该语句，若分支变量的值取值2，则执行指令组2，且结束该语句，······若分支变量不取所列出的值，则执行指令组k。
例：关于switch 和input 的用法。[useswitch.m]
⑸其它常用指令
●input指令提示用户从键盘输入数值、字符串或表达式，并接收该输入，语法为：
user_entry=input(‘message’)
user_entry=input(‘message’,’s’)
●pause 指令使程序运行暂停，语法为
pause:暂停执行程序，等待用户按任意键继续。
pause(n):使程序暂停n秒后继续执行。
例：for n=1:4
x=-1:0.1:1;y=x.^n;plot(x,y)
pause
end
●return 指令结束return 指令所在函数的执行，返回到主调函数或者命令窗口。
●break 指令中断执行，用在循环语句内表示跳出循环。
●error（’message’）显示出错信息，终止程序执行。
例：编写用区间迭代法求函数零点的程序。[intfzero.m]
4、其它一些有用的函数
●fcnchk 函数验证函数
f=fcnchk(fun),fun可以是由字符串表示的函数表达式，（这时返回一个inline函数），也可以是函数句柄，或是函数名字符串。
f=fcnchk(fun，‘vectorized’),生成向量化函数，例如用.*代替* 举例将intfzero.m 文件中的语句f=inline(fun);换成
f=fcnchk(fun);
●nargin ,nargout 函数中输入参数或输出参数的个数。
●% 注释语句
●find 寻找数组中非零元素对应的下标。
S=find(A),[I,J]=find(A)
例：重编上面的分段函数[myfun3.m]
●取整函数round(x),ceil(x),floor(x),fix(x)
第二部分Matlab作图
1、曲线图
●plot(x,y) 作出以数据(x(i),y(i))为节点的折线图，其中x,y为同维
数的向量。
●plot(x,y,s) 其中s是由颜色、标记、线型参数组成的字符串
颜色标记线型
b blue . point - solid
g green o circle : dotted
r red x x-mark -. dashdot
c cyan + plus -- dashed
m magenta * star (none) no line
y yellow s square
k black d diamond
w white v triangle (down)
^ triangle (up)
< triangle (left)
> triangle (right)
p pentagram
h hexagram
● plot(x1,y1,s1,x2,y2.s2,…) 在同一个坐标系中作出由向量对
(x1,y1),(x2,y2),…为节点的折线。
例：在同一坐标系中作出函数31y x x =--和0.2
sin(5)
y x
x =
在区间[-1,2]
上的图形。
● hold on (off) 保持（释放）图形窗口
● polar(theta,rho) 作以（theta,rho ）为坐标的极坐标图形，theta,rho
为同维数的向量
例：作出四叶玫瑰线4sin 2ρ
θ
=的图形。
● fplot(fun,[a,b]) 作出函数fun 在区间[a,b]上的图形。 ● plot3(x,y,z) 作空间曲线的图形，x,y ,z 为同维向量。 例：作出曲线sin ,cos ,x t t y t t z t ===的图形。
● subplot(m,n,k) 将图形窗口分成m n ?个子图形窗口，将当前操作
定位在第k 个子图形窗口。
2、曲面图
[X,Y]=meshgrid(x,y) 生成以数组x,y 为坐标的网格矩阵 mesh(X,Y ,Z) 绘制网面图，X 、Y 、Z 是同维矩阵 surf(X,Y ,Z) 绘制曲面图，与mesh 用法类似。 例：作出曲面2
2
x
y
z xe --=在22,22x y -≤≤-≤≤上的图形
● ezmesh(fun) 轻松绘出二元函数fun 的曲面图（easy to use mesh ） 绘图区域为[2,2,2,2]ππππ--
● ezmesh(fun,[xmin,xmax,ymin,ymax]) 在指定区域绘图 ● ezmesh(fun,…,’circ ’) 绘图区域为圆域
上机练习
1、设x 为一个长度为n 的数组，编程求下列均值和标准差
1
1,1n
i i x x s n n ===
>∑ 2、求满足0
ln(1)100m
n n =+>∑的最小m 值。
3、用循环语句形成Fibonacci 数列12
121,,3,4,k k k F F F F F k --===+=
，并
验证极限
1
12k k F F -+→。（提示：计算至两边误差小于精度810-）
4、分别用for 和while
循环结构编写程序，求出6
101i K ==∑。并考虑
一种避免循环语句的程序设计，比较不同算法的运行时间。 5、求出所有的“水仙花数”。所谓“水仙花数”是指一个三位数，其各位数字的立方和等于该数本身。例如，153是一个“水仙花数”，因为333153153=++。
6、假定某天的气温变化记录如下表，试作图描述这一天的气温变化规律。
7、作出下列函数图像
⑴曲线22
sin(2),22
y x x x x =---≤≤ （要求分别使用plot 或fplot 完
成） ⑵椭圆
2
2
14
9
x
y
+
=
⑶抛物面22,3,3z x y x y =+<<
⑷曲面422232226,
3,313z x x y x y x y x y =++---+<-<<
⑸空间曲线sin ,cos ,cos(2),02x t y t z t t π===<< ⑹半球面2sin cos ,2sin sin ,2cos ,02,0/2
x y z φθφθφθπφπ===≤≤≤≤
⑺三条曲线合成图123sin ,sin sin(10),sin ,0y x y x x y x x π
===-≤≤
8、作下列分段函数图
1.1 1.1 1.11.1 1.1
x y x x x ?
?>=≤-<- 9、用MA TLAB 函数表示下列函数，并作图
2222220.75 3.75 1.56
0.75 3.75 1.50.5457 1(,)0.7575e
1 1 0.5457e 1 y x x
y x y x x
e
x y p x y x y x y -------+??????
???
+>=-<+≤+≤-
MATLAB程序设计 用MATLAB语言编写的程序，称为M文件。M文件可以根据调用方式的不同分为两类：命令文件(Script File)和函数文件(Function File)。 例3-1 分别建立命令文件和函数文件，将华氏温度f转换为摄氏温度c。 程序1：首先建立命令文件并以文件名f2c.m存盘。 clear; %清除工作空间中的变量 f=input('Input Fahrenheit temperature：'); c=5*(f-32)/9 然后在MATLAB的命令窗口中输入f2c，将会执行该命令文件，执行情况为： Input Fahrenheit temperature：73 c =22.7778 程序2：首先建立函数文件f2c.m。 function c=f2c(f) c=5*(f-32)/9 然后在MATLAB的命令窗口调用该函数文件。 clear; y=input('Input Fahrenheit temperature：'); x=f2c(y) 输出情况为： Input Fahrenheit temperature：70 c =21.1111 x =21.1111 3.1.2 M文件的建立与打开 M文件是一个文本文件，它可以用任何编辑程序来建立和编辑，而一般常用且最为方便的是使用MATLAB提供的文本编辑器。
1．建立新的M文件 为建立新的M文件，启动MATLAB文本编辑器有3种方法： (1) 菜单操作。从MATLAB主窗口的File菜单中选择New菜单项，再选择M-file命令，屏幕上将出现MATLAB 文本编辑器窗口。 (2) 命令操作。在MATLAB命令窗口输入命令edit，启动MATLAB文本编辑器后，输入M文件的内容并存盘。 (3) 命令按钮操作。单击MATLAB主窗口工具栏上的New M-File命令按钮，启动MATLAB文本编辑器后，输入M文件的内容并存盘。 2．打开已有的M文件 打开已有的M文件，也有3种方法： (1) 菜单操作。从MATLAB主窗口的File菜单中选择Open命令，则屏幕出现Open对话框，在Open对话框中选中所需打开的M文件。在文档窗口可以对打开的M文件进行编辑修改，编辑完成后，将M文件存盘。 (2) 命令操作。在MATLAB命令窗口输入命令：edit 文件名，则打开指定的M文件。 (3) 命令按钮操作。单击MATLAB主窗口工具栏上的Open File命令按钮，再从弹出的对话框中选择所需打开的M文件。 3.2 程序控制结构 3.2.1 顺序结构 1．数据的输入 从键盘输入数据，则可以使用input函数来进行，该函数的调用格式为： A=input(提示信息，选项)； 其中提示信息为一个字符串，用于提示用户输入什么样的数据。 如果在input函数调用时采用's'选项，则允许用户输入一个字符串。例如，想输入一个人的姓名，可采用命令： xm=input('What''s your name?','s'); 2．数据的输出 MATLAB提供的命令窗口输出函数主要有disp函数，其调用格式为
矩量法m ａtla ｂ程序设计实例： Ha ｌlen 方程求对称振子天线 一、条件与计算目标 已知: 对称振子天线长为Ｌ,半径为a ，且天线长度与波长得关系为,，设,半径a=0、０００0001，因此波数为。 目标： 用H ａll ｅn 方程算出半波振子、全波振子以及不同值得对应参数值。 求:（1）电流分布 （2）E 面方向图 (二维）,H 面方向图（二维）,半波振子空间方向性图(三维） 二、对称振子放置图 图1 半波振子得电流 分布 半波振子天线平行于z 轴放置，在ｘ轴与ｙ轴上得分量都为零，坐标选取方式有两种形式，一般选取图1得空间放置方 式。图１给出了天线得电流分布情况，由图可知，当天线很细时,电流分布近似正弦分布。 三、Ha ｌｌen 方程 得解题思路 ()()()()2 1 ' ' ' ' 12,cos sin sin 'z z i z z z z i z k z G z z dz c kz c kz E k z z dz j ωμ'++=-?? 对于中心馈电得偶极子，Hallen 方程为 ()22'1222 ('),'cos sin sin ,2L L i L L V i z G z z dz c kz c kz k z z j η + -- ++= <<+? 脉冲函数展开与点选配，得到 ()1121 ,''cos sin sin ,1,2,,2n n N z i n m m m m z n V I G z z dz c kz c kz k z m N j η +''=++= =???∑? 上式可以写成 矩阵形式为 四、结果与分析 （1）电流分布
1.编写程序：计算1/3+2/5+3/7+……+10/21 法一： s=0; for i=1:10 s=s+i/(2*i+1); end s s = 4.4096 法二： sum((1:10)./(3:2:21)) ans = 4.4096 2.编写程序：计算1~100中即能被3整除，又能被7整除的所有数之和。 s=0; for i=1:100 if mod(i,3)==0&&mod(i,7)==0 s=s+i; end,end s s = 210 3.画出y=n！的图（1<=n<=10）,阶乘的函数自己编写，禁用MATLAB自带的阶乘函数。 x=1:10; for i=1:10 try y(i)=y(i-1)*i; catch y(i)=1; end,end plot(x,y)
12345678910 0.511.522.533.54x 10 6 4.一个数恰好等于它的因子之和，这个数就称为完数。例如，6的因子为1，2,3，而6=1+2+3，因此6就是一个完数。编程找出2000以内的所有完数。 g=[]; for n=2:2000 s=0; for r=1:n-1 if mod(n,r)==0 s=s+r; end end if s==n g=[g n]; end end g g =6 28 496
5.编写一个函数，模拟numel函数的功能，函数中调用size函数。 function y=numelnumel(x) m=size(x); y=m(1)*m(2); numelnumel([1 2 3;4 5 6]) ans = 6 6. 编写一个函数，模拟length函数的功能，函数中调用size函数。 function y=lengthlength(x) m=size(x); y=max(m(1),m(2)); lengthlength([1 2 3;4 5 6]) ans = 3 7.求矩阵rand（5）的所有元素和及各行平均值，各列平均值。 s=rand(5); sum=sum(sum(s)) mean2=mean(s,2) mean1=mean(s) sum = 13.8469
MATLAB 程序设计方法及若干程序实例 樊双喜 (河南大学数学与 信息科学学院开封475004) 摘要本文通过对 MATLAB 程序设计中的若干典型问题做简要的分析和总结,并在此基础上着重讨论了有关算法设计、程序的调试与测试、算法与程序的优化以及循环控制等方面的问题.还通过对一些程序实例做具体解析,来方便读者进行编程训练并掌握一些有关MATLAB 程序设计方面的基本概念、基本方法以及某些问题的处理技巧等.此外,在文章的最后还给出了几个常用数学方法的算法程序, 供读者参考使用.希望能对初学者进行 MATLAB 编程训练提供一些可供参考的材料,并起到一定的指导和激励作用,进而为MATLAB 编程入门打下好的基础. 关键字算法设计；程序调试与测试；程序优化；循环控制 1 算法与程序 1.1 算法与程序的关系算法被称为程序的灵魂,因此在介绍程序之前应先了 解什么是算法.所谓算 法就是对特定问题求解步骤的一种描述.对于一个较复杂的计算或是数据处理的问题,通常是先设计出在理论上可行的算法,即程序的操作步骤,然后再按照算法逐步翻译成相应的程序语言,即计算机可识别的语言. 所谓程序设计,就是使用在计算机上可执行的程序代码来有效的描述用于解决特定问题算法的过程.简单来说,程序就是指令的集合.结构化程序设计由于采用了模块分化与功能分解,自顶向下,即分而治之的方法,因而可将一个较复杂的问题分解为若干子问题,逐步求精.算法是操作的过程,而程序结构和程序流程则是算法的具体体现. 1.2MATLAB 语言的特点 MATLAB 语言简洁紧凑,使用方便灵活,库函数极其丰富,其语法规则与科技人员的思维和书写习惯相近,便于操作.MATLAB 程序书写形式自由,利用其丰富
matlab程序设计例题及答案 1.编写程序：计算1/3+2/5+3/7+……+10/21 法一： s=0; for i=1:10 s=s+i/(2*i+1); end ss = 法二： sum((1:10)./(3:2:21)) ans = 2.编写程序：计算1~100中即能被3整除，又能被7整除的所有数之和。 s=0; for i=1:100 if mod(i,3)==0&&mod(i,7)==0 s=s+i; end,end ss = 210 3.画出y=n！的图,阶乘的函数自己编写，禁用MATLAB 自带的阶乘函数。 x=1:10; for i=1:10 try y(i)=y(i-1)*i; catch y(i)=1; end,end plot(x,y) 10612345678910 4.一个数恰好等于它的因子之和，这个数就称为完数。
例如，6的因子为1，2,3，而6=1+2+3，因此6就是一个完数。编程找出20XX以内的所有完数。 g=; for n=2:20XX s=0; for r=1:n-1 if mod(n,r)==0 s=s+r; end end if s==n g=[g n]; end end g g =6 28 496 5.编写一个函数，模拟numel函数的功能，函数中调用size函数。 function y=numelnumel(x) m=size(x); y=m(1)*m(2); numelnumel([1 2 3;4 5 6]) ans = 6 6. 编写一个函数，模拟length函数的功能，函数中调用size函数。 function y=lengthlength(x) m=size(x); y=max(m(1),m(2)); lengthlength([1 2 3;4 5 6]) ans = 3