当前位置：文档库 › 第五节Matlab编程

第五节Matlab编程

§5.1关系运算§5.1.1 比较运算

§5.1.2逻辑运算

§5.2控制流§5.2.1 for 循环

§5.2.2 while 循环

§5.2.3 if和break语句§5.3 M文件、命令文件及函数文件§5.3.1 M 文件

§5.3.3函数文件

§5.3.2命令文件§5.4 字符串、输入及输出§5.4.1 echo、input、pause、keyboard§5.4.2串和宏串

§5.4.3外部程序

§5.4.4输入输出数据

§5.1关系运算

§5.1.1 比较运算

比较两个同阶矩阵有下面六种相关操作符：

相关操作符

< 小于

<=小于等于

> 大于

>=大于等于

==等于

~=不等于

表5.1.1.1相关操作符

比较两个元素的大小，结果是“1”表明为真，结果是“0”表明为假．

例如2+2~=4

结果是“0”，表明为假．

例如一个6阶魔术方阵，矩阵元素计算满足各种条件：

A=magic(6)

ans =

35 1 6 26 19 24

3 32 7 21 23 25

31 9 2 22 27 20

8 28 33 17 10 15

30 5 34 12 14 16

4 36 29 13 18 11

阶数为n的魔术方阵，即n×n矩阵，是由1～n2的整数组成(n=6)．仔细观察这个矩阵，我们会发现任何行和、任何列和都相等．另外，每个3×3子行列式的对角线元素和，都可被3整除．为了显示这一特性，键入：

p=(rem(A,3)==0)

p =

0 0 1 0 0 1

1 0 0 1 0 0

0 1 0 0 1 0

0 0 1 0 0 1

1 0 0 1 0 0

0 1 0 0 1 0

为了再仔细地观察这个模式，可以用format+格式画出矩阵的压缩格式．此格式用“+”代表正元素，“-”代表负元素，空格代表0．

format +

p =

+ +

find 函数在关系运算中很有用，它可以在0-1矩阵中找非零元素的下标．

若y是一个向量，例如：y=[1 3 2 4 3.5 2.9]，则find(y<3.0)，将指出y的分量在哪些位置上小于3.0．

ans = 1 3 6

即：向量y的第1、3、6位置上的元素小于3.0．

当输入x==NaN时结果为NaN，因为按照IEEE算法规定任何具有NaN的操作，结果都是NaN．调试NaN很有用，例如测试x，输入isnan(x)函数，如果x元素是不定值则得1，否则得0．isfinite(x)更有用，如-∞

§5.1.2逻辑运算

&与

|或

~非

表5.1.2.1逻辑运算符

“＆”和“|”操作符可比较两个标量或两个同阶矩阵．对于矩阵来说必须符合规则，如

果A和B都是0-1矩阵，则A＆B或A|B也都是0-1矩阵，这个0-1矩阵的元素是A和B对应元素之间逻辑运算的结果，逻辑操作符认定任何非零元素都为真，给出“1”，任何零元素都为假，给出“0”．

非(或逻辑非)是一元操作符，即~A：当A是非零时结果为“0”；当A为“0”时，结果为“1”．因此下列两种表示：

p | (~p) 结果为1．

p ＆(~p) 结果为0．

any和all函数在连接操作时很有用，设x是0-1向量，如果x中任意有一元素非零时，any(x)返回“1”，否则返回“0”；all(x)函数当x的所有元素非零时，返回“1”，否则也返回“0”．这些函数在if语句中经常被用到．如：

if all(A<5)

do something

end

§5.2控制流

Matlab与其它计算机语言一样，也有控制流语句．控制流语句可使原本简单地在命令行中运行的一系列命令或函数，组合成为一个整体——程序，从而提高工作效率．

§5.2.1 for 循环

Matlab与其它计算机语言一样有do或for循环，完成一个语句或一组语句在一定时间内反复运行的功能．例如：

for i = 1:n , x( i )=0, end

x的第一个元素赋0值，如果n＜1，结构上合法，但内部语句不运行，如果x不存在或比n元素小，额外的空间将会自动分配．

多重循环写成锯齿形是为了增加可读性．例如：

m=9;n=9;

for i = 1:m

for j=1:n

A( i, j ) = 1/( i + j - 1);

end

程序的说明：

（1）事实上，上述程序给出了Hilbert矩阵的构造过程，可参见函数hilb(n)．

（2）语句内部使用分号，表示计算过程不输出中间结果．

（3）循环后的A命令表示显示矩阵A的结果．

（4）每个for语句必须以end语句结束，否则是错误的．

for循环的通用形式为：

for v=expression

statements

end

其中expression表达式是一个矩阵，因为Matlab中都是矩阵，矩阵的列被一个接一个的赋值到变量v，然后statements语句运行．

通常expression是一些m:n或m:k:n仅有一行的矩阵，并且它的列是个简单的标量．但如注意到expression可以为矩阵，即v可以为向量，对某些问题的处理将大大简化．

§5.2.2 while 循环

Matlab中的while 循环语句为一个语句或一组语句在一个逻辑条件的控制下重复未知的次数．

它的一般形式为：

while expression

statements

end

当expression的所有运算为非零值时，statements语句组将被执行．如果判断条件是向量或矩阵的话，可能需要all或any函数作为判断条件．

例如计算expm(A)，在A并不是太大时，直接计算expm(A)是可行的．

expm(A)=I+A+A^2/2!+A^3/3!+…注意：这里的I表示单位矩阵．程序为：

E = 0*A;

F = E + eye(size(E)); N = 1;

while norm(F,1) > 0,

F = A*F/N;

E = E + F;

N = N + 1;

end

§5.2.3 if和break语句

下面介绍if语句的二个例子．

(1) 一个计算如何被分成三个部分，用符号校验：

if n<0

A=negative(n)

elseif mod(n,2)==0

A=even(n)

else

A=odd(n)

end

其中的三个函数negative(n)、even(n)、odd(n)是自编的输出函数．参见下面的函数文件．

(2) 这个例子涉及数论中一个很有趣的问题，取任何的正整数，如果是偶数，用2除；如果是奇数，用3乘，并加上1，反复这个过程，直到你的整数成为1．这个极有趣不可解的问题是：有使这个过程不中止的整数吗?

%classic "3n+1"problem from number theory

while 1

n=input('Enter n, negative quits: ');

if n<=0 break,end

while n>1

if rem(n,2) == 0 %% 是连续的2个等号

n=n/2

else

n=3*n+1

end;

end

这个过程能永远进行吗?

程序的说明：

（1）本程序用到了if语句与while语句，过程比较复杂；

（2）使用input函数，可使程序在执行过程中，从键盘输入一个数（矩阵）；

（3）break语句提供了程序跳出死循环的途径．

§5.3 M文件、命令文件及函数文件

§5.3.1 M文件

Matlab通常使用命令驱动方式，当单行命令输入时，Matlab立即处理并显示结果，同时将运行说明或命令存入文件．

Matlab语句的磁盘文件称作M文件，因为这些文件名的未尾是.M形式，例如一个文件名为bessel.m，提供bessel函数语句．

一个M文件包含一系列的Matlab语句，一个M文件可以循环地调用它自己．

M文件有两种类型：

第一类型的M文件称为命令文件，它是一系列命令、语句的简单组合．

第二类型的M文件称为函数文件，它提供了Matlab的外部函数．用户为解决一个特定问题而编写的大量的外部函数可放在Matlab工具箱中，这样的一组外部函数形成一个专用的软件包．

这两种形式的M文件，无论是命令文件，还是函数文件，都是普通的ASCII文本文件，可选择编辑或字处理文件来建立．

§5.3.2命令文件

当一个命令文件被调用时，Matlab运行文件中出现的命令而不是交互地等待键盘输入，命令文件的语句在工作空间中运算全局数据，对于进行分析解决问题及做设计中所需的一长串繁杂的命令和解释是很有用的．

例如：一个自编的命令文件fibo.m，用于计算Fibonnaci数列

% An M-file to calculate Fibonnaci numbers

f=[1, 1 ]; i = 1;

while f(i)+f(i+1)<1000

f(i+2)=f(i)+f(i+1);

i=i+1;

end

plot(f)

在Matlab命令窗口中键入fibo命令，并回车执行，将计算出所有小于1000的Fibonnaci 数，并绘出图形．

要注意的是：文件执行后，f和i变量仍然留在工作空间．

§5.3.3函数文件

如果M文件的第一行包含function，这个文件就是函数文件，它与命令文件不同，所定义变量和运算都在文件内部，而不在工作空间．函数被调用完毕后，所定义变量和运算将全部释放．函数文件对扩展Matlab函数非常有用．

例如：一个自编的函数文件mean.m，用于求向量的（或矩阵按列的）平均值function y=mean(x)

% MEAN Average or mean value，For Vectors,

% MEAN (x) returns the mean value

% For matrix MEAN (x) is a row vector

% containing the mean value of each column

[m,n]=size(x);

if m==1

m=n;

end

y=sum(x)/m;

磁盘文件中定义的新函数称为mean函数，它与Matlab函数一样使用，例如z为从1到99的实数向量：

z=1:99;

计算均值：mean(z)

ans=

mean.m程序的说明：

（1）第一行的内容：函数名，输入变量，输出变量，没有这行这个文件就是命令文件，而不是函数文件．

（2）％：表明％右边的行是说明性的内容注释．前一小部分行来确定M文件的注释，并在键入help mean后显示出来．显示内容为连续的若干个％右边的文字．（3）变量m,n和y是mean的局部变量，在mean运行结束后，它们将不在工作空间z 中存在．如果在调用函数之前有同名变量，先前存在的变量及其当前值将不会改变．

再例如：一个计算标准差的函数文件stat.m

function [mean,stdev]=stat(x)

[m,n]=size(x);

if m==1

m=n

end

mean=sum(x)/m;

stdev=sqrt(sum(x.^2)/m-mean.^2);

stat表明返回多输出变量是可能的．

又如：使用多输入变量计算矩阵秩函数

function r=rank(x,tol)

% rank of a matrix

s=svd(x);

if(nargin==1)

tol=max(size(x))*s(1)*eps;

end

r=sum(s>tol);

这个变量说明利用永久变量nargin确定输入变量的个数，变量nargout虽然这里没有使用，但它包含有输出变量的个数．

一些有用的说明：

当M函数文件第一次在Matlab运行时，它被编译并放入内存，以后使用时不用重新编译即可得到．

what命令：显示磁盘当前目录中的M文件，

dir命令：列出所有文件．

一般而言，输入一个名字到Matlab，例如键入whoopie命令，Matlab用以下步骤解释：

(1) 看whoopie是否为变量．

(2) 检验whoopie是否为在线函数．

(3) 检验whoopie文件的当前目录．

(4) 将whoopie看成Matlab的PATH中的一个文件，在Matlab PATH目录中搜索．

如果whoopie存在，Matlab首先将其作为变量而不是作为函数．

§5.4 字符串、输入及输出

§5.4.1 echo、input、pause、keyboard

一般来说，当一个M文件运行时，文件的命令不在屏幕上显示，而echo命令则使M文件运行时，命令在屏幕上显示，这对于调试、演示相当有用．

input功能：输入Input('How many apples')给用户一个提示串，等待，然后显示用户通过键盘输入的大量表达式．可以用input命令建立驱动M文件的菜单．

与input功能相同，但功能更强的keyboard命令将计算机作为一个命令文件来调用，放入M文件中，此特性对调试或正在运行期间修改变量很有用．

pause命令：使用户暂停运行一个程序，当再按任一键时恢复执行，pause(n)等待n秒钟后再继续执行．

§5.4.2串和宏串

字符串用单个引号输入到Matlab中，例如：

s='Hello'

结果显示为：

s =

Hello

字符存在向量中，每个元素就是一个字符，如：

size(s)

ans =

1 5

表明S为一个1×5的矩阵，有五个元素．字符以ASCII值存入，abs函数或double函数将显示以下值(即Hello的ASCII值)

abs(s)

ans =

72 101 108 108 111

getstr函数，使向量作为字符显示，而不显示ASCII值．

disp可在变量中显示字符．sprintf, num2str和int2str可以将数字转换成串．

字符变量通过括号连成大串．例如：

s='hello';

s=[s,' world']

s =

hello world

eval是与字符变量—起工作的函数，执行简单字符宏调用．eval( t )执行包含在t内的字符．如果t是任何Matlab表达式或语句的源字符，则字符串被解释执行．例如：t='eye(2)', eval(t)

结果为：

ans=

1 0

0 1

又例如，给矩阵元素赋值

t='1/(i+j-1)';

for i=1:n

for j=1:n

a( i, j)=eval(t);

end

这儿有一个例子，介绍如何一起使用eval与load命令，装入十个具有顺序文件名的文件中的数据：

fname='mydata';

for i = 1:10

eval([ 'load ', fname, int2str( i )])

end

§5.4.3外部程序

Matlab与外部独立程序的通讯方式可以是多种多样的，下面介绍其中的一个办法：

(1) Matlab中将变量存入磁盘

(2) 运行外部程序（读数据文件，进行处理），将结果写到磁盘上

(3) 将处理后的文件装回到工作空间中

例如：用外部程序gareqn找garfield方程的结果：

function y=garfield(a,b,q,r)

save gardata a,b,q,r

! gareqn

load gardata

使用FORTRAN或其它语言写gareqn程序，使其可以读gardata.mat，进行处理，将结果存入文件中．

这个程序可将计算机的“连接码”提供给Matlab，在许多系统中它将新的目标码连接到程序中比物理联接要方便得多．

§5.4.4输入输出数据

可使用各种方法将其它程序和外部世界的数据送入Matlab，同样可把Matlab数据输送

到外部世界，使你的程序以Matlab使用的文件形式直接计算数据．

最好的方法取决于多少数据，数据是否可读，什么形式等：

(1) 清晰的元素表输入：

如果你有少量数据，比如说小于10~15个元素，使用方括号[]输入．

(2) 使用文本编辑建立命令文件，将数据列为清晰的元素表输入．如果数据不是可读形式，又不得不以一种方法键入，可以重复运行M文件，重复修改数据．

(3) 如果数据以ASCII形式存贮，并有固定长度，行尾有回车符，各数间有空格的文件称为flat file(ASCII的flat file可由普通文本编辑来编辑)，flat file通过load命令直接读进Matlab，结果存入名为文件名的变量中去．

(4) 将数据文件译成Matlab文件形式，使用load命令，translate程序由Matlab中的应用程序库支持，translate程序将ASCII文件、二进制文件、FORTRAN非格式文件和DIF文件转换为Matlab使用的特定的MA T文件，当磁盘文件中存有大量数据时，这个方法输入最好．

Matlab数据输出到外部世界的方法：

(1) 小矩阵时：使用diary命令建立日志文件，在文件中列出变量，用文本编辑处理日志文件，日志的输出包括运行中的Matlab命令．

(2) 使用save命令存入变量，退出Matlab，用translate程序将MAT文件转换成任一种其它文件形式．

第六节Matlab符号运算

§6.1 Matlab符号运算的工作流程

§6.1.1工作过程

§6.1.2 核心工具

§6.1.3 符号变量确定原则

§6.2 Matlab的六大常见符号运算§6.2.1 因式分解§6.2.2 计算极限

§6.2.3 计算导数

§6.2.4 计算不定积分、定积分、反常积分

§6.2.5 符号求和§6.2.6 解代数方程和常微分方程

Matlab本身并没有符号计算功能，1993年通过购买Maple的使用权后，开始具备符号运算的功能．符号运算的类型很多，几乎涉及数学的所有分支．

§6.1 Matlab符号运算的工作流程

§6.1.1工作过程

§6.1.2 核心工具

sym函数VS syms 语句

sym函数：构造符号变量和表达式: a=sym('a')

（Construct symbolic numbers, variables and objects）

syms语句：构造符号对象的简捷方式（Short-cut for constructing symbolic objects）

§6.1.3 符号变量确定原则

（1）除了i 和j之外，字母位置最接近x的字母；若距离相等，则取ASCII码大的；

（2）若没有除了i 与j以外的字母，则视x为默认的符号变量；

（3）可利用函数findsym(string,N)来询问在众多符号中，哪N个为符号变量．例如：键入findsym(3*a*b+y^2,1)，即可得到答案y．更多的例子见下表：

符号表达式默认符号变量

a*x^2+b*x+c x

1/(4+cos(t))t

4*x/y x

2*a+b b

2*i x

§6.2 Matlab的六大常见符号运算

§6.2.1 因式分解

syms x

f=x^6+1;

s=factor(f)

结果为：

s=(x^2+1)*(x^4-x^2+1)

§6.2.2 计算极限

求极限：

syms h n x

L=limit('(log(x+h)-log(x))/h',h,0) %%单引号可省略掉M=limit('(1-x/n)^n',n,inf)

结果为：

L =1/x

M =exp(-x)

§6.2.3 计算导数

．syms a x; y=sin(a*x);

A=diff(y,x)

B=diff(y,a)

C=diff(y,x,2)

结果为：

A = cos(a*x)*a

B = cos(a*x)*x

C = -sin(a*x)*a^2

§6.2.4 计算不定积分、定积分、反常积分

，

．

f=(x^2+1)/(x^2-2*x+2)^2;

g=cos(x)/(sin(x)+cos(x));

h=exp(-x^2);

I=int(f)

J=int(g,0,pi/2)

K=int(h,0,inf)

结果为：

I =1/4*(2*x-6)/(x^2-2*x+2)+3/2*atan(x-1)

J =1/4*pi

K =1/2*pi^(1/2)

§6.2.5 符号求和

求级数的和S, 以及前十项的部分和S1．

syms n

S=symsum(1/n^2, 1, inf)

S1=symsum(1/n^2,1,10)

结果为：

S =1/6*pi^2

S1 =1968329/1270080

重要说明：当求函数项级数的和S2时，可用命令：

syms n x

S2=symsum(x/n^2, n, 1, inf)

S2 =1/6*x*pi^2

两点说明：

（1）注意观察S2与S1的细微区别！

（2）当通项公式的Matlab表达式较长时，表达式要加上单引号．后面的练习中会遇到此问题．

§6.2.6 解代数方程和常微分方程

利用符号表达式解代数方程所需要的函数为solve(f)，即解符号方程式f．

例如：求一元二次方程a*x^2+b*x+c=0的根．

f=sym('a*x^2+b*x+c') 或f='a*x^2+b*x+c'

solve(f)

ans=

[1/2/a*(-b+(b^2-4*c*a)^(1/2))]

[1/2/a*(-b-(b^2-4*c*a)^(1/2))]

solve(f, a)

ans=

利用符号表达式可求解微分方程的解析解，所需要的函数为dsolve(f)，使用格式：dsolve('equation1', ' equation2', …)

其中：equation为方程或条件．写方程或条件时，用Dy 表示y 关于自变量的一阶导数，用D2y 表示y 关于自变量的二阶导数，依此类推．

1．求微分方程的通解．

syms x y %定义x,y为符号

dsolve('Dy=x', 'x')

ans =

1/2*x^2+C1

试比较：

若写成：

syms x y %定义x,y为符号

dsolve('Dy=x')

结果将是什么？是否正确？为什么？

2．求微分方程的特解．

syms x y

dsolve('D2y=x+Dy', 'y(0)=1', 'Dy(0)=0', 'x')

ans =

-1/2*x^2+exp(x)-x

试比较：

若写成：

syms x y

dsolve('D2y=x+Dy', 'y(0)=1', 'Dy(0)=0')

结果将是什么？是否正确？为什么？

3．求微分方程组的通解．

syms x y

[x,y]=dsolve('Dx=y+x, Dy=2*x')

x =

1/3*C1*exp(-t)+2/3*C1*exp(2*t)+1/3*C2*exp(2*t)-1/3*C2*exp(-t)

y =

2/3*C1*exp(2*t)-2/3*C1*exp(-t)+2/3*C2*exp(-t)+1/3*C2*exp(2*t)

试比较：

若写成：

（1）dsolve('Dx=y+x, Dy=2*x')

结果将是：

ans =

x: [1x1 sym]

y: [1x1 sym]

试解释此结果的含义．

若写成：

（2）[x,y]=dsolve('Dx=y+x, Dy=2x')

结果将是：

x =

exp(t)*C1+C2*exp(t)-C2-2-2*t

y =

C2+2*t

是否正确？为什么？

实验一特殊函数与图形

一、问题背景与实验目的

二、相关函数（命令）及简介

三、实验内容

四、自己动手

一、问题背景与实验目的

著名的Riemann函数大家都很熟悉了，但是关于它的图像你是否清楚呢？除了最上面那几点，其他都很难画吧？你想不想看看下面那些“挤在一起”的点是怎样分布的呢？还有几何中的马鞍面、单叶双曲面等是怎样由直线生成的，是不是也想目睹一下呢？这些，都离不开绘图．

实际上绘图一直是数学中的一种重要手段，借助图形，往往可以化繁为简，

使抽象的对象得到明白直观的体现．比如函数的基本性质，一个图形常可以使之一目了然，非常有效．它虽不能代替严格的分析与证明，但在问题的研究过程中，可以帮助研究人员节约相当一部分精力．此外，它还可以使计算、证明、建模等的结果得到更明白易懂的表现，有时，这比科学论证更有说服力．同时，数学的教学与学习过程也离不开绘图．借助直观的图形，常可以使初学者更容易接受新知识．如数学分析中有不少函数，其解析式着实让人望而生畏，即使对其性质作了详尽的分析，还是感到难明就里；但如果能看到它的图形，再配合理论分析，则问题可以迎刃而解．又如在几何的学习中，会遇到大量的曲线与曲面，也离不开图形的配合．

传统的手工作图，往往费力耗时，效果也不尽理想．计算机恰恰弥补了这个不足，使你可以方便地指定各种视角、比例、明暗，从各个角度进行观察．本实验通过对函数的图形表示和几个曲面（线）图形的介绍，一方面展示它们的特点，另一方面，也将就Matlab软件的作图功能作一个简单介绍．大家将会看到，Matlab 的作图功能非常强大．

二、相关函数（命令）及简介

1．平面作图函数：plot，其基本调用形式：

plot(x,y,s)

以x作为横坐标，y作为纵坐标．s是图形显示属性的设置选项．例如：x=-pi:pi/10:pi;

y=sin(x);

plot(x,y,'--rh','linewidth',2,'markeredgecolor','b','markerfacecolor','g')

图1

在使用函数plot时，应当注意到当两个输入量同为向量时，向量x与y必须维数相同，而且必须同是行向量或者同是列向量．

绘图时，可以制定标记的颜色和大小，也可以用图形属性制定其他线条特征，这些属性包括：

linewidth 指定线条的粗细．

markeredgecolor 指定标记的边缘色

markerfacecolor 指定标记表面的颜色．

markersize 指定标记的大小．

若在一个坐标系中画几个函数，则plot的调用格式如下：

plot(x1,y1,s1,x2,y2,s2,……)

2．空间曲线作图函数：plot3，它与plot相比，只是多了一个维数而已．其调用格式如下：

plot3(x,y,z,s)．例如：

x=0:pi/30:20*pi;

y=sin(x);z=cos(x);

plot3(x,y,z)

得到三维螺旋线：

图2

3．空间曲面作图函数：

（1）mesh函数．绘制彩色网格面图形．调用格式：

mesh(z)，mesh(x,y,z)和mesh(x,y,z,c)．

其中，mesh(x,y,z,c)画出颜色由c指定的三维网格图．若x、y均为向量，则length(x)=n，length(y)=m，[m,n]=size(z)．

（2）surf在矩形区域内显示三维带阴影曲面图．调用格式与mesh类似．（3）ezmesh用符号函数作三维曲面网格图．

调用格式：

ezmesh(x,y,z)

其中x = x(s,t), y = y(s,t),z = z(s,t)．画图区域默认为：-2*pi < s < 2*pi 且-2*pi < t < 2*pi.

或者用格式：

ezmesh(x,y,z,[smin,smax,tmin,tmax])

（4）ezsurf用符号函数作三维曲面图．调用格式与ezmesh类似．

（5）sphere画球体命令．

4．meshgrid，调用格式：

[x,y]=meshgrid(m,n)，

这里的m，n为给定的向量，可以定义网格划分区域和划分方法．矩阵x和矩阵y是网格划分后的数据矩阵．

5．图像的修饰与其他函数：

（1）axis equal 控制各个坐标轴的分度，使其相等；

（2）colormap设置绘图颜色．

调用格式：

colormap([r g b])

其中r,g,b都是0-1之间的数．

或者用格式：

colormap(s)

s为颜色映像．下面举几个常用的例子：

颜色映像相应的颜色系颜色映像相应的颜色系

autumn 红黄色系hsv 色调饱和色系

gray 线性灰色系hot 黑红黄白色系

cool 青和洋红色系pink 柔和色系

（3）grid网格函数grid on添加网格．grid off取消网格．

（4）find找出符合条件的元素在数组中的位置．调用格式：

y=find(条件)

例如：输入：

a=[4 5 78 121 4 665 225 4 1];

b=find(a>7)

输出： b =3 4 6 7

三、实验内容

数学分析中，特别是积分部分，我们接触了不少有趣的函数，由于其中有的不是一一对应的，用上面的方法无法画出它们的图像，这时就只能用参数了．此外还有些图形只能用参数来画，比如空间曲线，在计算机上不接受“两个曲面的交线”这种表示，所以也只能用参数来实现．

用参数方式作图的关键在于找出合适的参数表示，尤其是不能有奇点，最好也不要用到开方．所以要找的参数最好是有几何意义的．当然这也不可一概而论，需要多积累经验．

1．利用函数plot在一个坐标系中画以下几个函数图像，要求采用不同颜色、不同线形、不同的符号标记．函数为：．

程序如下：

t=0:pi/20:2*pi;

x=sin(t);

y=cos(t);

z=sin(2*t);

plot(t, x, '--k*', t, y, '-rs', t, z, ':bo')

图像如下：

图3

2．绘制类似田螺线的一条三维螺线(方程自己设计)．

程序如下：

t=0:.1:30;

x=2*(cos(t)+t.*sin(t));

y=2*(sin(t)-t.*cos(t));

z=1.5*t;

plot3(x,y,-z) %取–z 主要是为了画图看起来更清楚

axis equal

图像如下：

图4

3．利用函数，绘制一个墨西哥帽子的图形．程序如下：

[a,b]=meshgrid(-8:.5:8); %先生成一个网格

c=sqrt(a.^2+b.^2)+eps;

z=sin(c)./c;

mesh(a,b,z)

axis square

图像如下：

图5

思考：这里的eps 是什么？其作用是什么？

4．利用surf绘制马鞍面图形（函数为：）．

程序如下：

[x,y]=meshgrid(-25:1:25,-25:1:25);

z=x.^2/9-y.^2/4;

surf(x,y,z)

title('马鞍面')

grid off

图像如下：

图6

5．分别用ezmesh和ezsurf各绘制一个圆环面，尝试将两个圆环面放在一个图形

界面内，观察它们有什么不同之处．

提示：圆环面的方程为：，而圆环面的参数方程为：

程序参见附录1．

图像如下：

图7

6．绘制黎曼函数图形，加深对黎曼函数的理解．

说明：黎曼函数的定义为

程序参见附录2．

图像如下：