计算机建模与仿真-第3章.ppt

第三章,MATLAB数值运算,标识符：标识符是标志变量名、常量名、函数名和文件名的字符串的总称。MATLAB中的标识符最长允许31个字符，合法字符是52个英文字母（大小写看成不同字符）、10个数字和下划线。第一个字符必须是英文字母。,3.1 MATLAB数值运算基础,例:MYvar12,My_var12,MyVar12_均为有效的变量名 _MyVar12为无效的变量名 ABC和ABc表达的是不同的变量,变量及其赋值数,MATLAB中的实数只有一种数据格式，那就是双精度格式（即64位二进制或8个字节），它的有效值是十进制16位，动态范围是10的308次幂。数的存储和运算全按照同一格式进行，这样使得编程简单，减少差错，改善了人机交互，但使得运行速度和内存空间的消耗上付出了代价。此为MATLAB的一大特色。,变量及其赋值矩阵及其元素的赋值,格式1：变量=表达式（或常量）格式2：表达式例:a=15+202-50+3*9a=392 b=30 b=30 a/b%表达式不赋给变量时，MATLAB自动给一个临时标量ans来暂时存储。ans=13.0667,这里表达式可进行+、*、/、sqrt、log、exp,变量及其赋值矩阵及其元素的赋值,分号和逗号的区别 x=2,y=3%逗号隔开，屏幕有回显x=2y=3 m=2;n=3;%分号隔开，无回显 m%在提示符后直接输入变量名可查看变量的值m=2,变量及其赋值矩阵及其元素的赋值,矩阵有如下运算：,注意:4/2 42含义一样么？,变量及其赋值矩阵及其元素的赋值,通过直接输入矩阵的元素构造矩阵：用中括号 把所有矩阵元素括起来同一行的不同数据元素之间用空格或逗号间隔用分号（；）指定一行结束可分成几行进行输入，用回车符代替分号数据元素可以是表达式，系统将自动计算结果,变量及其赋值矩阵及其元素的赋值,例：输入矩阵A、B的值 A=1 2 3 4;5 6 7 8;9 10 11 12;13 14 15 16 B=1,sqrt(25),9,132,6,10,7*23+sin(pi),7,11,154,abs(-8),12,16,变量及其赋值矩阵及其元素的赋值,变量中的元素所在位置用（）中的数字来注明，以逗号隔开，用户可单独对其赋值。如a（3，2）=8，表示a矩阵中的第三行第二列元素值为8。若赋值元素的下标超出了原来矩阵的大小，矩阵的行列会自动扩展。修改矩阵A中元素的数值 A=1 2 3 4;5 6 7 8;9 10 11 12;13 14 15 16;A(1,1)=0;A(2,2)=A(1,2)+A(2,1);A(4,4)=cos(0);A(5,5)=8;A=0 2 3 4 0 5 7 7 8 0 9 10 11 12 0 13 14 15 1 0 0 0 0 0 8,变量及其赋值矩阵及其元素的赋值,给全行元素赋值可用“：”需要抽去某行元素另构成一个矩阵用“”A=1 2 3 4;5 6 7 8;9 10 11 12;13 14 15 16;A(4,:)=3 4 5 6 A(1,2,:)=%亦可表示为A(1:2,:)=A=9 10 11 12 3 4 5 6,变量及其赋值矩阵及其元素的赋值,注意：空矩阵就是指没有元素的矩阵。对任何一个矩阵赋值，就是使它的元素都消失掉。这完全不同于“零矩阵”，后者是元素存在，只是其数值为零而已。空矩阵是使矩阵缩短时不可缺少的矩阵。,变量及其赋值矩阵及其元素的赋值,矩阵另外的赋值方式：矩阵下标与子矩阵提取A(m,n)提取第m行，第n列元素A(:,n)提取第n列元素A(m,:)提取第m行元素A(m1:m2,n1:n2)提取第m1行到第m2行和第n1列到 第n2列的所有元素A(m1,m2,n1,n2)提取第m1行和第m2行，第n1列和 第n2列的所有元素A(m:end,n)提取从第m行到最末行和第n列的子块A(:)得到一个长列矢量，该矢量的元素按 矩阵的列进行排列,变量及其赋值矩阵及其元素的赋值,a=1 2 3 4;5 6 7 8;9 10 11 12;13 14 15 16;b=a(3,4)%提取第3行，第4列元素 c=a(:,1)%提取第1列元素 d=a(1,:)%提取第1行元素 e=a(1:3,1:3)%提取第1行到第3行和第1列到第3列的所有元素 f=a(1,3,1,3)%提取第1行和第3行，第1列和第3列的所有元素 g=a(2:end,3)%提取从第2行到最末行和第3列的子块 h=a(:)%得到一个长列向量，该向量的元素按矩阵的列进行排列 h=a(:)%得到一个长行向量,虚数符号MATLAB启动时定为i、j，在写复数的虚部时可以不用乘号，连写在数字后面。如3+5i等。如果用户在程序中另外给i、j赋值，则他们的虚数意义就失效。此时需用clear i,j将原来赋给i，j的值清除。复数赋值格式：格式1：按元素逐个赋值格式2：将其实部和虚部分别赋值,变量及其赋值复数,c=5+2i z=1+2i,3+4i;5+6i,7+8i z=1,3;5,7+2,4;6,8*i%此处*不可以省,变量及其赋值复数,MATLAB中所有的运算符和函数都对复数有效。conj(X)表示为共轭，即把X的虚部反号。如Z=3+5i，conj（Z）=3-5i。是共轭转置运算符，对实矩阵把行号与列号交换，对复矩阵除行列交换外，还要把矩阵元素取共轭。若求复数的转置而不要共轭，增将conj和结合起来使用。,变量及其赋值复数,f=sqrt(1+2i)g=sin(3+5i)f*f,f*g z=1+2i,3+4i;5+6i,7+8i w=z,u=conj(z),v=conj(z),变量及其赋值复数,变量检查命令：who，whos clear%删除工作区中所有定义过的变量 whos%查看当前工作区内变量信息，无显示表示没有定义的变量 xy=1;yx=2;%对变量赋值 whos Name Size Bytes Class xy 1x1 8 double array yx 1x1 8 double arrayGrand total is 2 elements using 16 bytes whoYour variables are:xy yx clear xy yx%删除变量xy及yx whos xy%这时变量xy已经不存在了?Undefined function or variable xy.,变量及其赋值变量检查,变量及其赋值变量检查,在检查中不显示内部变量(具体见P12页表2.1中的特殊变量和函数）,变量及其赋值变量检查,注意：在其他语言中若遇到inf 或nan这些特殊变量时，系统就停止运算退出。而MATLAB不停止运算，仍给结果赋予inf 或nan，并继续把程序执行完。这是MATLAB的又一大好处，可避免因为一个数据不好而破坏全局。出现inf 或nan后，系统对他做任何运算，结果仍为inf 或nan，此为IEEE标准。,变量及其赋值基本赋值矩阵,变量及其赋值基本赋值矩阵,f1=ones(3,2)f2=zeros(2,3)f3=eye(3)f4=rand(3,5),f5=randn(3,5)f6=linspace(0,1,5)%在0和1中均匀的生成5个点（包括0和1），形成5维向量。f7=logspace(0,1,5)%在100和101中按对数均匀的生成5个点（包括0和1），形成5维向量。,变量及其赋值矩阵的组合与显示,大矩阵可由若干个小矩阵组成，但必须其行列数正确，恰好填满全部元素。一个矩阵中所有元素用同一显示格式。有一个是小数则都显示为小数。当矩阵中的最大元素小于 0.001，或其最小元素大于1000时，MATLAB 会把其公因子提出来。元素大小差别很大时，不宜用矩阵显示。,变量及其赋值矩阵的组合与显示,f1=ones(3,2),f2=zeros(2,3),f3=magic(3),f4=eye(2)f5=linspace(0,1,5)fb1=f1,f3;f4,f2 fb2=fb1;f5 f=0.000073 5.33e-6,3.2 MATLAB的矩阵运算,（1）矩阵的加减乘法（2）矩阵的除法和线形方程解（3）矩阵的乘方和幂次运算（4）矩阵结构形式的提取与变换,两个矩阵进行加减的条件是其行数和列数对应相同，简称阶次相同 size(A)=size(B)矩阵的加减即其对应元素进行加减。如果其中一个变量是标量，则会自动扩展成各元素均取此标量值的同阶次矩阵。,3.2 矩阵的初等运算矩阵的加减法,x=1 2 3;4 5 6;7 8 9;y=9 8 7;6 5 4;3 2 1;mx,nx=size(x)my,ny=size(y)x+y x-y x-1%1是标量，则会自动扩展成各元素均取此标量值的33矩阵,3.2 矩阵的初等运算矩阵的加减法,3.2 矩阵的初等运算矩阵的加减法,注意：size和lengthsize有两个输出量。格式：m,n=size(x)，其中m和n分别表示矩阵x的行数和列数。length只有一个输出量，格式：n=length(x)，n表示矩阵x的列数，或行向量x的列数例：x=1 2 3;4 5 6 mx,nx=size(x)m=length(x)y=1 2 3 4 5 6 7 8 n=length(y),3.2 矩阵的初等运算矩阵的乘法,两矩阵相乘C=A*B的条件是A的列数等于B的行数，简称内阶数相同。如矩阵A为n*p阶，矩阵B为p*m阶，则C是n*m阶。内阶数不同的矩阵不能相乘。矩阵的乘法不符合交换律：A*B B*A，左乘不等于右乘。若一个变量是标量，则按标量乘法进行运算。,3.2 矩阵的初等运算矩阵的乘法,A=1,2,3;4,5,6,B=2,4,0;1,3,5,C=1,4,7;8,5,2,;3,6,0,D=1,2,3;4,5,6;7,8,9 A*B%内阶数不同，无法进行运算 A*B A*B C*D D*C%上两句表示乘法不符合交换律 C*2%若一个变量是标量，则按标量乘法进行运算,3.2 矩阵的初等运算矩阵的乘法,线性方程组D*X=B，如果D非奇异，即它的逆矩阵inv（D）存在；则其解用MATLAB表示为：X=inv(D)*B=DB符号称为左除，即分母放在左边。左除的条件：B的行数等于D的行数。A=1,2,3;4,5,6,D=1,4,7;8,5,2,;3,6,0 DA,DA,3.2 矩阵的初等运算矩阵的除法,3.2 矩阵的初等运算矩阵的除法,线性方程组X*D1=B1，如果D1非奇异，即它的逆矩阵inv（D1）存在；则其解用MATLAB表示为：X=B1*inv(D1)=D1/B1符号/称为右除。右除的条件：B1的列数等于D1的列数。A=1,2,3;4,5,6,D=1,4,7;8,5,2,;3,6,0 A/D,A/D,3.2 矩阵的初等运算线性方程组的求解,3.2 矩阵的初等运算线性方程组的求解,MATLAB的乘幂函数“”、指数函数expm、对数函数logm、开方函数sqrtm是对矩阵进行运算，相当于矩阵整体运算。另外，矩阵的乘法和除法运算也是把矩阵作为一个整体来运算。,3.2 矩阵的初等运算矩阵的乘方和幂次运算,s=1 2;3 4,D=1 4 7;8 5 2;3 6 0 u1=sqrtm(s)v1=expm(s)logm(D)D2,3.2 矩阵的初等运算矩阵的乘方和幂次运算,除此之外，MATLAB还可对矩阵中的元素分别进行运算。此称为元素群运算。如“.”，“.*”，“.”，exp，log，sqrt等。（2.3节介绍）,3.2 矩阵的初等运算矩阵的乘方和幂次运算,在做矩阵运算时，往往需要提取其中的某些特殊结构的元素来组成新的矩阵，有时则需要改变矩阵的排列。除了2.1节讲过的提取行、列和将行列转置的语句外，MATLAB还提供了些改变矩阵结构的函数。如fliplr，flipud，reshape，rot90，diag，tril，triu等。,3.2 矩阵的初等运算矩阵结构形式的提取与变换,A=8 1 6 0;3 5 7 1;4 9 2 2 B=fliplr(A)B=flipud(A)B=reshape(A,2,6)B=rot90(A)B=diag(A)B=tril(A)B=triu(A)B=A(:),3.2 矩阵的初等运算矩阵结构形式的提取与变换,3.2 MATLAB的数组运算,（1）数组及其赋值（2）数组运算的四则运算和幂次运算（3）数组函数 数组运算能大大简化编程，提高运算的效率，是MATLAB优于其他语言的一大特色。,数组通常是指单行或单列的矩阵，也称为向量。其赋值方法：（1）逐个输入 a=1 3 9 10 15 16%采用空格和逗号分隔构成行向量 b=1;3;9;10;15;16%采用分号隔开构成列向量（2）x=初值：增量：终值 增量为1时，则赋值可改为x=初值：终值（3）线性分割函数x=linspace（起点，终点，点数）（4）对数分割函数x=logspace（起点对数，终点对数，点数）,3.2 数组运算数组及其赋值,a=1 3 9 10 15 16 b=1;3;9;10;15;16 t=0:0.02:1 z=10:-3:-5 k=1:6 theta=linspace(0,2*pi,9)w=logspace(0,1,11),3.2 数组运算数组及其赋值,数组运算：就是把两矩阵按逐个元素进行运算（不按照线性代数的规则）。为了与矩阵作为整体的运算符号相区别，要在运算符“*、/、”前加一点符号“.”。另外还有exp，log，sqrt，用于进行元素群运算数组运算的英文原文是（ArrayAlgorithm），,3.2 数组运算数组的四则和幂次运算,X=1,2,3,Y=4,5,6 Z=X.*Y Z=X.Y Z=X./Y Z=X.Y Z=X.2 2.X Y s=1 2;3 4,D=1 4 7;8 5 2;3 6 0 u2=sqrt(s)v2=exp(s)log(D),3.2 数组运算数组的四则和幂次运算,3.2 数组运算数组的四则和幂次运算,元素群和矩阵幂次的比较 D=1,4,7;8,5,2,;3,6,0 Z=D3 D.3 Z=3.D,D=1 4 7 8 5 2 3 6 0Z=627 636 510 804 957 516 486 612 441ans=1 64 343 512 125 8 27 216 0Z=3 81 2187 6561 243 9 27 729 1,大多数MATLAB函数都适用于作数组运算，只有专门说明的几个除外，就是*、/、运算符和P51页的矩阵函数，具体见P49-表3-7、3-8,3.3 数组函数、矩阵函数,利用元素群造表 x=0:0.05:pi/4;x,sin(x),cos(x),tan(x)disp(x sin(x)cos(x)tan(x);disp(x,sin(x),cos(x),tan(x),3.2数组运算数组的函数运算,3.4 多项式函数库（polyfun）,多项式是形如P(x)=a0 xn+a1xn-1+an-1x+an的式子。在MATLAB中，多项式用行向量表示：P=a0 a1 an-1 an如：x4-34x3+5x2-1多项式可表示为1-34 5 0-1,3.4 多项式函数库（polyfun）,多项式的行向量除了直接输入外，还可由poly（a）命令构成，根据a的形式不同，有两种处理结果：若a是方阵，poly（a）生成的多项式是矩阵a的特征多项式若a是形如 r0 r1 rn-1 rn的向量，命令poly（a）会生成多项式(x-r0)(x-r1)(x-rn-1)(x-rn)所对应的n阶多项式a0 xn+a1xn-1+an-1x+an。的行向量 a0 a1 an-1 an,a=2 9 0 0;0 4 1 4;7 5 5 1;7 8 7 4 pa=poly(a)ppa=poly2str(pa,x)%（向量显示成多项式形式）b=4 6 8 9 pb=poly(b)ppb=poly2str(pb,y),3.4 多项式函数库（polyfun）,3.4 多项式函数库（polyfun）,多项式的四则运算多项式求导、求根和求值多项式拟合多项式插值线性微分方程的解（residue）,3.4 多项式函数库（polyfun）多项式的四则运算,多项式的四则运算参与加减运算的多项式应该具有相同的阶次，若阶数不同，低阶的多项式必须用零填补至高阶多项式的阶次。多项式乘法采用conv函数除法由deconv函数完成。,a=2,4,6,8 b=3,6,9 c=a+0,b d=a-0,b e=conv(a,b)q,r=deconv(c,a)%q为商,r为余数,3.4 多项式函数库（polyfun）多项式的四则运算,3.4 多项式函数库（polyfun）多项式求导、求根和求值,多项式求导、求根和求值polyder（求导）roots（求根）polyval（求值）,a=2,4,6,8 b=3,6,9 p=polyval(b,1)%当x=1时的值 p=polyval(b,2,4)%当x=2、4时的值 x=1 1;1 1;p=polyval(b,x)%当x为某个方阵时该多项式的值 q=roots(a)pa=poly2str(a,x)e=polyder(a)pa1=poly2str(e,x),3.4 多项式函数库（polyfun）多项式求导、求根和求值,3.4 多项式函数库（polyfun）多项式拟合,多项式拟合在实际工程应用中，我们往往测得一些分散的数据点，为了找出他们的内在规律性，就需要对这些分散的数据点运用最小二乘法、多项式或其他已知函数等方法来生成一个新的多项式或函数来逼近已知点，此为拟合。MATLAB具有强大的计算功能和绘图功能，所以可方便的进行曲线拟合并绘制拟合曲线。,多项式拟合曲线拟合涉及两个问题：最佳拟合是什么意思？应该用什么样的曲线进行拟合？最佳拟合被解释为：在数据点的最小误差平方和。且所用的曲线限定为多项式时，此时曲线拟合是相当简捷的。多项式的最小二乘曲线拟合。,3.4 多项式函数库（polyfun）多项式拟合,多项式拟合p=polyfit(x,y,n)x和y为已知N个数据点坐标向量。n是用来拟合的多项式的次数。p是求出的多项式的系数。n次多项式应该有n+1个系数，故p的长度为n+1。采用最小二乘法拟合。,3.4 多项式函数库（polyfun）多项式拟合,x=0:0.1:1;y=-.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2;subplot(2,3,1);plot(x,y,o)xi=linspace(0,1);a1=polyfit(x,y,1);y1=polyval(a1,xi);subplot(2,3,2);plot(x,y,o,xi,y1,b)a2=polyfit(x,y,2);y2=polyval(a2,xi);subplot(2,3,3),plot(x,y,o,xi,y2,r)a3=polyfit(x,y,3);y3=polyval(a3,xi);subplot(2,3,4),plot(x,y,o,xi,y3,r)a9=polyfit(x,y,9);y9=polyval(a9,xi);subplot(2,3,5),plot(x,y,o,xi,y9,b)a10=polyfit(x,y,10);y10=polyval(a10,xi);subplot(2,3,6),plot(x,y,o,xi,y10,b),3.4 多项式函数库（polyfun）多项式拟合,结果:a1=10.3185 1.4400a2=-9.8108 20.1293-0.0317a3=16.0758-33.9245 29.3246-0.6104a9=1.0e+005*-0.2531 0.9746-1.4494 0.9933-0.2387-0.0680 0.0501-0.0095 0.0008-0.0000a10=1.0e+006*-0.4644 2.2965-4.8773 5.8233-4.2948 2.0211-0.6032 0.1090-0.0106 0.0004-0.0000,3.4 多项式函数库（polyfun）多项式拟合,3.4 多项式函数库（polyfun）多项式插值,多项式插值其目的也是为了找到分散点的内在规律。而多项式插值时，其插值函数通常是分段的，人们不关心函数的表达式，而是关心插值得到的数据点。,多项式插值一维插值函数：yi=interp1(x,y,xi,methord)x和y为给定点，xi为插值范围内的任意点集的x坐标，yi是插值后的对应数据点集的y坐标。Methord为插值函数类型选项，有linear（一次线形，默认）,cubic（三次）,cubic spline（三次样条）等二维插值函数：zi=interp1(x,y,z,xi,yi,methord)例题略,3.4 多项式函数库（polyfun）多项式插值,3.4 多项式函数库（polyfun）线性微分方程的解（residue）,设微分方程为y+5y+4y+7y=3u+0.5u+4u求输入信号u为脉冲和阶跃信号时方程的解析解y用Laplace变换，可写出,因为脉冲输入u(s)=1阶跃输入u(s)=1/s所以y(s)可以表示成两个s的多项式之比,3.4 多项式函数库（polyfun）线性微分方程的解,3.4 多项式函数库（polyfun）线性微分方程的解,用r，p=residue(b，a)求出留数r和极点p它的反变换为：,a=1,5,4,7;b=3,0.5,4;%脉冲响应时参数 r,p=residue(b,a)t=0:0.2:10;yi=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t)+r(3)*exp(p(3)*t);subplot(1,2,1),h=plot(t,yi);a=1,5,4,7,0;b=3,0.5,4;%阶跃响应时参数 r,p=residue(b,a)ys=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t)+r(3)*exp(p(3)*t)+r(4);subplot(1,2,2),h=plot(t,ys);,3.4 多项式函数库（polyfun）线性微分方程的解,应用例题,(12+2*(7-4)/33 ans=0.6667,例：求解算术表达式的值,应用例题,a=5.67;b=7.811;exp(a+b)/log10(a+b)ans=6.3351e+005,例：计算下式的结果，其中a=5.67,b=7.811,应用举例,应用举例,作业,上机做习题集中第三章作业。,