matlab在科学计算中的应用02.ppt

第二章 MATLAB 语言程序设计基础MATLAB 语言的简洁高效性MATLAB 语言的科学运算功能MATLAB 语言的绘图功能MATLAB 庞大的工具箱与模块集MATLAB 强大的动态系统仿真功能 MATLAB MATLAB 语言是当前国际上自动控制领域的首选语言是当前国际上自动控制领域的首选计算机语言，也是很多理工科专业最适合的计算机数学计算机语言，也是很多理工科专业最适合的计算机数学语言。通过学习可更深入理解和掌握数学问题的求解思语言。通过学习可更深入理解和掌握数学问题的求解思想，提高求解数学问题的能力，为今后其他专业课程的想，提高求解数学问题的能力，为今后其他专业课程的学习提供帮助。学习提供帮助。MATLAB MATLAB语言的优势语言的优势:本章主要内容MATLAB 程序设计语言基础基本数学运算MATLAB语言流程控制MATLAB 函数的编写二维图形绘制三维图形绘制2.1 Matlab 的启动Windows Systems：u D双击桌面上的 MATLAB 图标 u 从“开始”中选择“MATLAB”启动程序打开Matlab的命令窗口。每个命令行提示符为.“help”的使用在命令行中键入 help matlab/general-General purpose commands.matlab/ops-Operators and special char.matlab/lang-Programming language const.matlab/elmat-Elementary matrices and ma.matlab/elfun-Elementary math functions.matlab/specfun-Specialized math functions.或者点击命令窗口中的“help选项”2.2 MATLAB 程序设计语言基础MATLAB 语言的变量命名规则命名规则是：（1）变量名必须是不含空格的单个词；（2）变量名区分大小写；（3）变量名最多不超过19个字符；（4）变量名必须以字母打头，之后可以是 任意字母、数字或下划线，变量名中 不允许使用标点符号例：NetCost,Left2Pay,x3,X3,z25c5 是允许的变量名；Net-Cost,2pay,%x,sign 是错误的变量名注意：不能使用Matlab保有的变量名，如pi，log等MATLAB 的保留常量数学运算符号及标点符号数学运算符号及标点符号Examples:2+3/4*5ans=5.7500Matlab 中公式计算顺序：1.quantities in brackets,2.powers 2+32(=2+9)=11,3.*/,working left to right(3*4/5=12/5),4.+-,working left to right(3+4-5=7-5),符号“e”用与表示很大或很小的数：-1.3412e+03=-1.3412 103-1341.2 -1.3412e-01=-1.3412 10(-1)=-0.13412 format 命令控制输出结果的长度。format long,format short(e)（1）MATLAB的每条命令后，若为逗号或无逗号或无标点标点符号，则显示命令的结果；若命令后为分分号号，则禁止显示结果.（2）“%”后面所有文字为注释.（3）“.”表示续行.双精度数值变量IEEE标准，64位(占8字节)，11指数位，53数值位和一个符号位 double()函数的转换其他数据类型uint8()，无符号8位整形数据类型，值域为0至255，常用于图像表示和处理。（节省存储空间，提高处理速度）int8(),int16(),int32(),uint16(),uint32()数值型数据结构符号型，sym(A),常用于公式推导、解析解解法 符号变量声明 syms var_list var_props 例：syms a b real(注意变量中间用空格分隔)syms c positive 符号型数值可采用变精度函数求值(variable-precision arithmetic)vpa(A),或 vap(A,n)vpa(pi)ans=3.1415926535897932384626433832795 vpa(pi,60)ans=3.14159265358979323846264338327950288419716939937510582097494符号型变量数据类型字符串型数据：用单引号括起来。多维数组：是矩阵的直接扩展，多个下标。单元数组：将不同类型数据集成到一个变量名下面，用表示；例：用Ai,j可表示单元数组A的第i行,第j列的内容。类与对象：允许用户自己编写包含各种复杂详细的变量，可以定义传递函数。MATLAB支持的其它数据结构直接赋值语句 赋值变量赋值表达式 例：a=pi2 a=9.8696例：表示矩阵 B=1+9i,2+8i,3+7j;4+6j 5+5i,6+4i;7+3i,8+2j 1iB=1.0000+9.0000i 2.0000+8.0000i 3.0000+7.0000i 4.0000+6.0000i 5.0000+5.0000i 6.0000+4.0000i 7.0000+3.0000i 8.0000+2.0000i 0+1.0000i B=1+9i,2+8i,3+7j;4+6j 5+5i,6+4i;7+3i,8+2j 1i;%不显示结果，但赋值MATLAB 的基本语句结构函数调用语句返回变量列表函数名（输入变量列表）例：a,b,c=my_fun(d,e,f,c)冒号表达式 v=s1:s2:s3 该函数生成一个行向量v，其中s1是起始值，s2是步长（若省略步长为1），s3是最大值。例：用不同的步距生成(0,p)间向量。v1=0:0.2:piv1=Columns 1 through 9 0 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 1.4000 1.6000 Columns 10 through 16 1.8000 2.0000 2.2000 2.4000 2.6000 2.8000 3.0000 v2=0:-0.1:pi%步距为负，不能生成向量，得出空矩阵v2=Empty matrix:1-by-0 v3=0:pi%默认步长为1v3=0 1 2 3 v4=pi:-1:0 逆序排列构成新向量v4=3.1416 2.1416 1.1416 0.1416 v5=0:0.4:pi,piv5=0 0.4000 0.8000 1.2000 1.6000 2.0000 2.4000 2.8000 3.1416基本语句格式 B=A(v1,v2)v1、v2分别表示提取行（列）号构成的向量。例：A=1,2,3,4;3,4,5,6;5,6,7,8;7,8,9,0A=1 2 3 4 3 4 5 6 5 6 7 8 7 8 9 0 B1=A(1:2:end,:)提取全部奇数行、所有列。B1=1 2 3 4 5 6 7 8子矩阵提取 B2=A(3,2,1,2,3,4)提取3，2，1行、2，3，4列构成子矩阵。A=B2=1 2 3 4 6 7 8 3 4 5 6 4 5 6 5 6 7 8 2 3 4 7 8 9 0 B3=A(:,end:-1:1)将A矩阵左右翻转，即最后一列排在最前面。B3=4 3 2 1 6 5 4 3 8 7 6 5 0 9 8 7矩阵表示矩阵转置数学表示(若A有复数元素，先转置再取各元素共轭复数值，Hermit转置）MATLAB 求解 BA C=A2.3 基本数学运算矩阵的代数运算矩阵加减法 C=A+B D=A-B注意维数是否相等注意其一为标量的情形矩阵乘法数学表示MATLAB 表示 C=A*B注意两个矩阵相容性 矩阵除法矩阵左除：AX=B，求 XMATLAB 求解：X=AB若A为非奇异方阵，则 X=A-1B最小二乘解（若A不是方阵）矩阵右除：XA=B，求 X MATLAB求解：X=B/A若A为非奇异方阵，则 X=BA-1最小二乘解（若A不是方阵）矩阵翻转左右翻转 B=fliplr(A)上下翻转 C=flipud(A)旋转 90o (逆时针）D=rot90(A)如何旋转180o？D=rot180(A)?Undefined function or variable rot180.D=rot90(rot90(A)矩阵乘方 A 为方阵，求 MATLAB 实现：F=Ax点运算-矩阵对应元素的直接运算数学表示:MATLAB 实现：C=A.*B例:A=1,2,3;4,5,6;7,8,0;B=A.A ()B=1 4 27 256 3125 46656 823543 16777216 1 C=A.*A ()C=1 4 9 16 25 36 49 64 0 a=1:5,b=6:10,a./ba=1 2 3 4 5b=6 7 8 9 10ans=Columns 1 through 3 0.16666666666667 0.28571428571429 0.37500000000000 Columns 4 through 5 0.44444444444444 0.50000000000000 a./aans=1 1 1 1 1逻辑变量：当前版本有逻辑变量对 double 变量来说，非 0 表示逻辑 1逻辑运算（相应元素间的运算）与运算 A&C或运算 A|C非运算 A异或运算 xor(A,C)矩阵的逻辑运算各种允许的比较关系 ,=,AA=1 2 3 4 5 6 7 8 0 find(A=5),大于或等于5元素的下标 ans=3 5 6 8矩阵的比较运算 i,j=find(A=5);i,j 显示行标，列标ans=A=3 1 1 2 3 2 2 4 5 6 3 2 7 8 0 2 3 all(A=5)某列元素全大于或等于5时，相应元素为1，否则为0。ans=0 0 0 any(A=5)某列元素中含有大于或等于5时，相应元素为1，否则为0。ans=1 1 1Size of Matrix A=5 7 9;1-3-7 ans=2 3 size(ans)ans=1 2 r,c=size(A)r=2 c=3Transpose of a matrix D=1 2 3 4 5;6 7 8 9 10;11 13 15 17 19;Dans=1 6 11 2 7 13 3 8 15 4 9 17 5 10 19 size(D),size(D)ans=3 5ans=5 3The Identity Matrix I=eye(3),x=8;-4;1,I*xI=1 0 0 0 1 0 0 0 1x=8 -4 1ans=8 -4 1Diagonal Matrices D=-3 0 0;0 4 0;0 0 2%维数比较高的时候，不现实D=-3 0 0 0 4 0 0 0 2 d=-3 4 2,D=diag(d)On the other hand,if A is any matrix,the command diag(A)extracts its diagonal entries:F=0 1 8 7;3-2-4 2;4 2 1 1,diag(F)F=0 1 8 7 3 -2 -4 2 4 2 1 1ans=0 -2 1Building Matrices(可以应用“小”的矩阵合成“大”的矩阵)C=0 1;3-2;4 2;x=8;-4;1;G=C xG=0 1 8 3 -2 -4 4 2 1 J=1:4;5:8;9:12;20 0 5 4;K=diag(1:4)J;J zeros(4,4)K=1 0 0 0 1 2 3 4 0 2 0 0 5 6 7 8 0 0 3 0 9 10 11 12 0 0 0 4 20 0 5 4 1 5 9 20 0 0 0 0 2 6 10 0 0 0 0 0 3 7 11 5 0 0 0 0 4 8 12 4 0 0 0 0解析结果的化简与变换MATLAB 实现：s1=simple(s)从各种方法中自动选择最简格式 s1,how=simple(s)化简并返回实际采用的化简方法 其中，s为原始表达式，s1为化简后表达式，how为采用的化简方法。其他常用化简函数（信息与格式可用 help命令得出）collect()合并同类项 expand()展开多项式 factor()因式分解 numden()提取多项式的分子和分母 sincos()三角函数的化简例：syms s;P=(s+3)2*(s2+3*s+2)*(s3+12*s2+48*s+64)P=(s+3)2*(s2+3*s+2)*(s3+12*s2+48*s+64)simple(P)%一系列化简尝试，得出计算机认为的最简形式ans=(s+3)2*(s+2)*(s+1)*(s+4)3 a,m=simple(P)%返回化简方法为因式分解方法，用 factor()函数将得同样结果 a=(s+3)2*(s+2)*(s+1)*(s+4)3m=factor expand(P)ans=s7+21*s6+185*s5+883*s4+2454*s3+3944*s2+3360*s+1152变量替换 其中，f为原表达式，用x*替换x得出新的。例：求其 Taylor 幂级数展开 syms a b c d t;%假设这些变量均为符号变量 f=cos(a*t+b)+sin(c*t)*sin(d*t);%定义给定函数 f(t)f1=subs(f,a,b,c,d,t,0.5*pi,pi,0.25*pi,0.125*pi,4)f1=-1.0000基本数论运算下取整、上取整、四舍五入、离0近方向取整、最简有理数、求模的余数、最大公约数、最小公倍数、质因数分解、判定是否为质数例：对下面的数据进行取整运算 -0.2765,0.5772,1.4597,2.1091,1.191,-1.6187 A=-0.2765,0.5772,1.4597,2.1091,1.191,-1.6187;floor(A)%向-inf 方向取整ans=-1 0 1 2 1 -2 ceil(A)%向+inf 方向取整ans=0 1 2 3 2 -1 round(A)%取最近的整数ans=0 1 1 2 1 -2 fix(A)%向 0 的方向取整ans=0 0 1 2 1 -1例：3x3 Hilbert 矩阵，试用 rat()函数变换 A=hilb(3);n,d=rat(A)将元素变换成最小有理数，n,d分别为分子、分母矩阵。n=1 1 1 1 1 1 1 1 1d=1 2 3 2 3 4 3 4 5例:1856120，1483720，最大公约数、最小公倍数，质因数分解。format long m=1856120;n=1483720;gcd(m,n),lcm(m,n)求m,n的最大公约数、最小公倍数。ans=1.0e+009*0.00000196000000 1.40508284000000 factor(lcm(n,m)对lcm(n,m)进行质因数分解。ans=2 2 2 5 7 7 757 947 n=sym(1483720);m=sym(1856120);gcd(m,n),lcm(m,n)ans=1960,1405082840 factor(lcm(m,n)ans=(2)3*(5)*(7)2*(757)*(947)例：1-100间质数 A=1:10;isprime(A)若向量A中某个整数值为质数，则相应位置为1，其他为零。ans=0 1 1 0 1 0 1 0 0 0 A=1:100;B=A(isprime(A)B=Columns 1 through 16 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 Columns 17 through 25 59 61 67 71 73 79 83 89 97rem(A,C)%A中元素对C中元素求模得出的余数。A=1:10;C=5;B=rem(A,C)B=1 2 3 4 0 1 2 3 4 0循环结构、条件结构、开关结构以及试探结构循环结构for 结构while 结构2.4 MATLAB 语言流程控制例：用循环求解 s=0;for i=1:100 s=s+i;end s=0;i=1;while(i sum(1:100)ans=5050例：用循环求解求最小的 m,s.t.s=0;m=0;while(s tic,s=0;for i=1:100000,s=s+1/2i+1/3i;end;tocelapsed_time=0.360000 tic,i=1:100000;s=sum(1./2.i+1./3.i);toc 向量化所需时间少elapsed_time=0.234000 i=1:10;s=1./2.i+1./3.i,ss=sum(1./2.i+1./3.i)s=0.8333 0.3611 0.1620 0.0748 0.0354 0.0170 0.0083 0.0041 0.0020 0.0010ss=1.4990转移结构例：用for循环求解求最小的 m，s.t.s=0;for i=1:10000 s=s+i;if s10000,break;endend ii=141开关结构和 C 语言的区别当开关表达式的值等于某表达式，执行该语句后结束该结构，不用 break当需要在开关表达式满足若干个表达式之一时执行某一程序段，则用单元形式（用大括号把这些表达式括起来，用逗号分隔）otherwise 语句，不是C语言中的 default（但与之等价）程序的执行结果和各个case顺序无关case 语句中条件不能重复，否则列在后面的条件将不能执行全新结构（首先试探性执行语句1，若执行过程中有错，将错误信息赋给保留的lasterr变量，并终止这段语句的执行，转而执行语句2。）应将不保险但快的算法放在语句1，保险的放在语句2；或语句2说明语句1失效原因。试探结构函数是 MATLAB 编程的主流方法除了函数外，还可以采用 M-script（M脚本文件）文件M-script 适合于小规模运算例：求满足 下面不等式的最小整数m.若最大值不为 10000，需修改程序对 m 和 10000 值的设置，不适合于M-script输入和数即可返回m值Matlab函数2.5 MATLAB 函数的编写 s=0;m=0;while(s=10000),m=m+1;s=s+m;end2.5.1 MATLAB 语言函数的基本结构nargin，nargout分别表示输入和返回变量的实际个数，此为MATLAB保留变量，只要进入该函数，MATLAB就将自动生成这两个变量。varargin,varargout 输入、输出变量列表（可变输入输出个数）。Matlab函数由function语句引导例：前面的要求，m,10000，满足：function m,s=findsum(k)s=0;m=0;while(s m1,s1=findsum(145323)m1=539s1=145530 无需修改程序例：编写过程中注意实现一下几点：若只给出一个输入参数，则会自动生成一个方阵在函数中给出合适的帮助信息检测输入和返回变量的个数edit myhilbfunction A=myhilb(n,m)%产生A=MYHILB(N,M)或A=MYHILB(N)；if nargout1,error(Too many output arguments.);endif nargin=1,m=n;elseif nargin=0|nargin2 error(Wrong number of input arguments.);End%对输入输出变量检测A1=zeros(n,m);%生成一个n行m列零矩阵for i=1:n for j=1:m A1(i,j)=1/(i+j-1);end,end 矩阵元素赋值if nargout=1,A=A1;elseif nargout=0,disp(A1);end help myhilb 产生A=MYHILB(N,M)或A=MYHILB(N)；A=myhilb(3,4)A=1.0000 0.5000 0.3333 0.2500 0.5000 0.3333 0.2500 0.2000 0.3333 0.2500 0.2000 0.1667 A=myhilb(4)A=1.0000 0.5000 0.3333 0.2500 0.5000 0.3333 0.2500 0.2000 0.3333 0.2500 0.2000 0.1667 0.2500 0.2000 0.1667 0.1429 A=myhilb(3,4,5)?Error using=myhilbToo many input arguments.例：函数的递归调用：阶乘function k=my_fact(n)if nargin=1,error(输入变量个数错误，只能有一个输入变量);endif nargout1,error(输出变量个数过多);endif abs(n-floor(n)eps|n1%如果 n1,进行递归调用 k=n*my_fact(n-1);elseif any(0 1=n)%0!=1!=1 k=1;end my_fact(11)ans=39916800其实MATLAB提供了求取阶乘的函数factorial(),其核心算法为 prod(1:n),从结构上更简单、直观，速度也更快。prod(1:11)%求1,2,11中各元素的乘积ans=39916800 prod(1:3:11)%数组1,4,7,10的乘积ans=280例：conv()可以计算两个多项式的积用 varargin 实现任意多个多项式的积function a=convs(varargin)a=1;for i=1:length(varargin),a=conv(a,varargini);end P=1 2 4 0 5;Q=1 2;F=1 2 3;D=convs(P,Q,F)D=1 6 19 36 45 44 35 30 poly2sym(D)ans=x7+6*x6+19*x5+36*x4+45*x3+44*x2+35*x+30 conv(Q,F)ans=1 4 7 62.5.2 可变输入输出个数Convolue：数学中卷积的运算 E=conv(conv(P,Q),F)%若采用 conv()函数，则需要嵌套调用E=1 6 19 36 45 44 35 30 poly2sym(E)ans=x7+6*x6+19*x5+36*x4+45*x3+44*x2+35*x+30 G=convs(P,Q,F,1,1,1,3,1,1)G=1 11 56 176 376 578 678 648 527 315 902.5.3 inline 函数和匿名函数inline 函数，可以免去文件f=inline(sin(x.2+y.2),x,y)MATLAB 7.0 inline一般用于描述某个数学函数方便例：函数2.6 二维图形绘制2.6.1 二维图形绘制基本语句构造向量:扩展扩展:绘制y矩阵每行和t关系的曲线(m条)条件:y矩阵的列数等于t向量的长度绘制 t 矩阵每行和 y 矩阵对应行关系的曲线条件:t 矩阵的行、列数与y矩阵的行、列数分别相等 t,y 都是矩阵例：选项为红色点划线且每个转折点用五角星表示r-.pentagram注:grid on(grid off),在图形上添加(取消)网格线hold on(hold off),保护(取消保护)当前坐标系title()，图形添加标题 xlabel()(ylabel(),给x(y)坐标轴添加标注例：x=-pi:0.05:pi;%以 0.05 为步距构造自变量向量 y=sin(tan(x)-tan(sin(x);%求出各个点上的函数值 plot(x,y)plot(x,y,r-.pentagram)x=-pi:0.05:-1.8,-1.801:.001:-1.2,-1.2:0.05:1.2,.1.201:0.001:1.8,1.81:0.05:pi;%以变步距方式构造自变量向量 y=sin(tan(x)-tan(sin(x);%求出各个点上的函数值 plot(x,y)%绘制曲线例：x=-2:0.02:2;%生成自变量向量 y=1.1*sign(x).*(abs(x)1.1)+x.*(abs(x)plot(x,y)plot(-2,-1.1,1.1,2,-1.1,-1.1,1.1,1.1)图形元素属性获取与修改 图形中，每条曲线、坐标轴、图形窗口分别是一个对象。可用set()函数设置对象的属性，用get()函数获得对象的某个属性。这两个语句在界面编程中特别有用。2.6.2 其他二维图形绘制语句 二维条形图、罗盘图、羽毛状图、直方图、极坐标图、阶梯图形、x-半对数图、彗星状轨迹图、误差限图形、二维填充图、对数图、磁力线图、火柴杆图、y-半对数图。例：绘制极坐标曲线 theta=0:0.01:6*pi;rho=5*sin(4*theta/3);polar(theta,rho)rho=5*sin(theta/3);polar(theta,rho)例：用不同曲线绘制函数表示正弦曲线 t=0:.2:2*pi;y=sin(t);%先生成绘图用数据 subplot(2,2,1),stairs(t,y)%分割窗口，在左上角绘制阶梯曲线 subplot(2,2,2),stem(t,y)%火柴杆曲线绘制 subplot(2,2,3),bar(t,y)%条型图绘制 subplot(2,2,4),semilogx(t,y)%横坐标为对数的曲线subplot(m,n,k)m:窗口分割的行数窗口分割的行数n:窗口分割的列数窗口分割的列数k:图形放置的位置图形放置的位置2.6.3 隐函数绘制及应用隐函数 ezplot(x2*sin(x+y2)+y2*exp(x+y)+5*cos(x2+y)x自选例：ezplot(x2*sin(x+y2)+y2*exp(x+y)+5*cos(x2+y),-10 10)2.7 三维图形绘制三维曲线绘制stem3（三维火柴杆型曲线）,fill3（三维填充图形）,bar3（三维直方图）等。例：参数方程 t=0:.1:2*pi;%构造 t 向量，注意下面的点运算 x=t.3.*sin(3*t).*exp(-t);y=t.3.*cos(3*t).*exp(-t);z=t.2;plot3(x,y,z),grid on%三维曲线绘制 stem3(x,y,z);hold on;plot3(x,y,z),grid off2.7.2 三维曲面绘制一般曲面绘制mesh()绘制网格图,surf()绘制表面图。其他函数，光照下 surfl(),等高线surfc()，瀑布型waterfall()等高线绘制 contour(),contour3()例:Butterworth 滤波器 x,y=meshgrid(0:31);n=2;D0=200;D=sqrt(x-16).2+(y-16).2);%求距离 z=1./(1+D.(2*n)/D0);mesh(x,y,z),%计算并绘制滤波器 axis(0,31,0,31,0,1)%重新设置坐标系，增大可读性 surf(x,y,z)%绘制三维表面图 contour3(x,y,z,30)三维等高线图，30等高线条数例：试绘制出二元函数 x,y=meshgrid(-2:.1:2);z=1./(sqrt(1-x).2+y.2)+1./(sqrt(1+x).2+y.2);Warning:Divide by zero.Warning:Divide by zero.(Type warning off MATLAB:divideByZero to suppress this warning.)surf(x,y,z),shading flat%修饰其显示形式 xx=-2:.1:-1.2,-1.1:0.02:-0.9,-0.8:0.1:0.8,0.9:0.02:1.1,1.2:0.1:2;yy=-1:0.1:-0.2,-0.1:0.02:0.1,0.2:.1:1;x,y=meshgrid(xx,yy);z=1./(sqrt(1-x).2+y.2)+1./(sqrt(1+x).2+y.2);surf(x,y,z),shading faceted;set(gca,zlim,0,20)%获得当前坐标轴对象的句柄例:Butterworth 滤波器 x,y=meshgrid(0:31);n=2;D0=200;D=sqrt(x-16).2+(y-16).2);z=1./(1+D.(2*n)/D0);%计算滤波器 subplot(221),surf(x,y,z),axis(0,31,0,31,0,1);view(0,90);%俯视图 subplot(222),surf(x,y,z),axis(0,31,0,31,0,1);view(90,0);%侧视图%subplot(223),surf(x,y,z),axis(0,31,0,31,0,1);view(0,0);%正视图 subplot(224),surf(x,y,z),axis(0,31,0,31,0,1);%三维图x-轴轴y-轴轴z-轴轴视角视角方位角方位角仰角仰角x-轴轴y-轴轴z-轴轴视角视角方位角方位角仰角仰角