《矩阵、数组和符号运算2new.ppt》由会员分享，可在线阅读，更多相关《矩阵、数组和符号运算2new.ppt（68页珍藏版）》请在三一办公上搜索。
1、,第3讲 矩阵、数组和符号运算2第3讲 Matlab流程控制结构3,第3讲矩阵、数组和符号运算,掌握内容：（1）了解 MATLAB 的符号变量，掌握 MATLAB 符号表达式、符号矩阵的两种创建方法。（2）掌握 MATLAB 符号数学函数的创建。（3）掌握符号矩阵的基本运算及MATLAB 关于不同精度的控制方法。（4）掌握符号微积分内容，包括求函数的极限、对符号表达式求导数和微分、符号积分、符号求和、傅立叶变换及其逆变换等。（5）掌握各种符号方程的求解方法和函数命令。（6）了解 MATLAB 可视化的符号函数分析界面及使用。（7）初步了解 MAPLE 的符号资源。,二、符号及运算,第3讲矩阵、
2、数组和符号运算,抽象运算：公式推导、因式分解、求解代数方程或微分方程的精确解 符号数学工具箱 1）通过基本符号数学工具箱的专用函数；符号表达式和符号矩阵的操作；多项式的化简、展开和代入；线性代数；微积分；符号方程的求解；特殊的数学函数。2）通过 maple.m、mpa.m 两个专门设计的 M 文件进行符号运算；3)通过 MATLAB 中的函数计算器（Function Caculator）。,第3讲矩阵、数组和符号运算,1、符号变量的创建a.符号矩阵定义函数 sym 函数 S=sym(arg)，从表达式 arg 创建一个 sym 对象 Sx=sym(x)x=sym(x,real)x=sym(x,
3、unreal)附加属性pi=sym(pi)delta=sym(1/10),第3讲矩阵、数组和符号运算,A=2/5,4/0.78,sqrt(23)/3;0.33,0.3333,log(4)输入数值矩阵AA=0.4000 5.1282 1.5986 0.3300 0.3333 1.3863 FA=sym(A)将数值矩阵A转化为符号矩阵FA FA=2/5,200/39,sqrt(23/9)不管数值矩阵的元素是以分数或是浮点数表示，转换后的符号矩阵都将以最接近有理式的形式给出。b.符号定义函数 syms 函数 syms arg1 arg2.syms a b c x y,2、符号表达式和矩阵的创建a.字
4、符串直接输入创建 符号表达式和符号方程对空格很敏感。因此，在创建符号表达式或符号方程时，不要在字符间任意加空格符；符号计算中出现的数字也是当作符号处理的；f=a*x2+b*x+cf=a*x2+b*x+c f=a*x2+b*x+c=0f=a*x2+b*x+c=0,第3讲矩阵、数组和符号运算,第3讲矩阵、数组和符号运算,这种方法输入符号矩阵与字符串矩阵的输入相似。但要保证在同一列中各元素字符串有同样的长度，在较短的字符串前后用空格符填充；这种方法要求符号矩阵每一行的两端都有方括号，而字符串矩阵仅在首尾有方括号。B=4+x x2 x;x3 5*x-3 x*aB=4+x x2 x x3 5*x-3 x
5、*a,第3讲矩阵、数组和符号运算,b.由 sym 命令创建 f=sym(a*x2+b*x+c)f=a*x2+b*x+c f1=sym(a*x2+b*x+c=0)f1=a*x2+b*x+c=0 A=sym(4+x,x2,x;x3,5*x-3,x*a)A=4+x,x2,x x3,5*x-3,x*a,第3讲矩阵、数组和符号运算,c.由 syms 命令创建 syms x a b c f=a*x2+b*x+cf=a*x2+b*x+c syms x a B=4+x x2 x;x3 5*x-3 x*aB=4+x,x2,x x3,5*x-3,x*a,不能创建符号方程 符号方程：equ1=(a*x+b=0),基
6、本运算 四则运算两个符号矩阵的大小相等方可进行加减运算，符号矩阵和符号标量的加减运算按照数组运算规则进行；两个符号矩阵只有内积相等时才可以进行乘法运算；符号的乘方运算 Sp，若 S 为符号表达式，p 可以为符号表达式或数值表达式；若 S 为符号矩阵，则 p 必须是整数。,第3讲矩阵、数组和符号运算,3、符号矩阵的运算,第3讲矩阵、数组和符号运算,a=sym(1/x,1/(x+1);1/(x+2),1/(x+3)a=1/x,1/(x+1)1/(x+2),1/(x+3)b=sym(x,1;x+2,0)b=x,1 x+2,0 b-aans=x-1/x,1-1/(x+1)x+2-1/(x+2),-1/
7、(x+3),abans=-6*x-2*x3-7*x2,3/2*x2+x+1/2*x3 6+2*x3+10*x2+14*x,-1/2*x3-2*x2-3/2*x a.bans=x2,x+1(x+2)2,0 a2ans=1/x2+1/(x+1)/(x+2),1/x/(x+1)+1/(x+1)/(x+3)1/(x+2)/x+1/(x+3)/(x+2),1/(x+1)/(x+2)+1/(x+3)2 exp(b)ans=exp(x),exp(1)exp(x+2),1,b符号简化Matlab提供了多种化简符号表达式的方法：降幂排列法（collect）展开法（expand）重叠法（horner）因式分解法（
8、factor）一般化简法（simplify）不定化简法（simple）分子分母形式法（numden）,第3讲矩阵、数组和符号运算,（1）降幂排列法collect（又称合并同类项）格式 R=collect(S)%对于多项式S中的每一函数，collect(S)按缺省变量x的次数合并系数。R=collect(S,v)%对指定的变量v计算，操作同上。,syms x y;R1=collect(exp(x)+x)*(x+2)R2=collect(x+y)*(x2+y2+1),y)R3=collect(x+1)*(y+1),x+y),R1=x2+(exp(x)+2)*x+2*exp(x)R2=y3+x*y2
9、+(x2+1)*y+x*(x2+1)R3=(y+1)*x+y+1,x+y,已知表达式，将f按x进行降幂排列，将g按 降幂排列,（2）展开法格式 R=expand(S)%对符号表达式S中每个因式的乘积进行展开计算。该命令通常用于计算多项式函数、三角函数、指数函数与对数函数等表达式的展开式。,E1=expand(x-2)*(x-4)*(y-t)E2=expand(cos(x+y)E3=expand(exp(a+b)3)E4=expand(log(a*b/sqrt(c),E1=x2*y-x2*t-6*x*y+6*x*t+8*y-8*tE2=cos(x)*cos(y)-sin(x)*sin(y)E3=
10、exp(a3)*exp(a2*b)3*exp(a*b2)3*exp(b3)E4=log(a*b/c(1/2),（3）重叠法格式 R=horner(P)%若P为一符号多项式的矩阵，该命令将矩阵的每一元素转换成嵌套形式的表达式R。,syms x yH1=horner(2*x4-6*x3+9*x2-6*x-4)H2=horner(x2+x*y;y3-2*y),H1=-4+(-6+(9+(-6+2*x)*x)*x)*xH2=x2+x*y；(-2+y2)*y,（4）因式分解法（factor）格式 factor(s)%符号表达式s的因式分解函数，说明：S为符号矩阵或符号表达式，常用于多项式的因式分解。,将
11、x 9-1分解因式syms x factor(x9-1),ans=(x-1)*(x2+x+1)*(x6+x3+1),（5）一般化简法（simplify）格式 R=simplify(S)说明：充分考虑符号表达式的各种运算、特殊函数的运算性质，经计算机比较后，给出认为表达式相对简单的一种化简方法。,syms x a b cR1=simplify(sin(x)4+cos(x)4)R2=simplify(exp(c*log(sqrt(a+b)S=(x2+5*x+6)/(x+2),sqrt(16);R3=simplify(S),R1=2*cos(x)4+1-2*cos(x)2R2=(a+b)(1/2*c
12、)R3=x+3,4,试编制程序，判断以下等式是否相等,（6）不定化简法（simple）simple(s)%s是矩阵或表达式，给出最简形式R,how=simple(s)%R为返回的最简形，how为简化过程中使用的主要方法。最笨拙，但综合前面5种化简方法的优点。用户可以根据所列出的结果进行比较和筛选。将表达式s的长度化到最短。若还想让表达式符合人的书写习惯，可使用函数Pretty。,不定化简法经常用的转化方法有：Combine(trig):以三角函数的运算性质对主对角代数式进行化简；Convert（exp）：将代数式尽量转化为由exp（x）、exp（ix）表示的指数形式；Convert（sinco
13、s）：将代数式尽量转化为由sin（x）、cos（x）表示的式子；Congvert（tan）：将代数式尽量转化为由tan（x）表示的式子,第3讲矩阵、数组和符号运算,因式分解 syms x factor(x9-1)ans=(x-1)*(x2+x+1)*(x6+x3+1)符号矩阵展开 syms x y expand(x+1)3)ans=x3+3*x2+3*x+1 expand(sin(x+y)ans=sin(x)*cos(y)+cos(x)*sin(y),同类项合并 syms x y collect(x2*y+y*x-x2-2*x,x)ans=(y-1)*x2+(y-2)*x分式通分 syms x
14、 y n,d=numden(x/y+y/x)n=x2+y2d=y*x,（7）分子分母法（numden）N,D=numden(A)将符号或数值矩阵A中的每一元素转换成整系数多项式的有理式形式，其中分子与分母是相对互素的。输出的参量N为分子的符号矩阵，输出的参量D为分母的符号矩阵。,syms x y a b c d;n1,d1=numden(sym(sin(4/5)n2,d2=numden(x/y+y/x)A=a,1/b;1/c d;n3,d3=numden(A),n2=x2+y2d2=y*xn3=a,1；1,dd3=1,b；c,1,C.符号变量（表达式）的提取与代入（1）subexpr 格式 Y
15、,SIGMA=subexpr(X,SIGMA)说明 找出符号表达式 X中相同的子表达式，将X中相同的、比较复杂的子字符串用符号SIGMA代替。,t=solve(a*x3+b*x2+c*x+d=0);r,s=subexpr(t,s);pretty(t)pretty(r)输入命令，并观察相同因子的提取,（2）符号变量的代入 格式subs(S,old,new)用符号new 代替表达式S 中的符号old Syms x ysubs(x+y,x,4)ans=4+ysyms a bf=subs(cos(a)+sin(b),a,b,sym(alpha),2)f=cos(alpha)+sin(2),试将以下表达
16、式中的变量x分别用y和mnt代替。f=a*x2+b*x+c,syms x y a b cf=a*x2+b*x+csf1=subs(f,x,y)sf2=subs(f,x,m+n*t),a.符号极限limit(F,x,a)：求当x-a时表达式F的极限值limit(F,a)，变量为由findsym定义的默认变量limit(F,x,a,right)/limit(F,x,a,left)（右/左极限值）Limit(F),a=0：求系统默认变量0时表达式F的极限值 syms x a t h limit(sin(x)/x)ans=1 limit(1+2*t/x)(3*x),x,inf)ans=exp(6*t)
17、,limit(1/x,x,0,right)ans=inf,6、符号微积分Matlab自变量确定原则：除i、j外，字母位置最接近x的小写字母为自变量；如果表达式中没有变量，x会被视为默认的变量。由函数findsym可一找到默认变量,练习：求下列函数极限,syms x hf=(sin(x+h)-sin(x)/hlimit(f,h,0),cos(x),b.符号积分int(S)不定积分int(S,v)int(S,a,b)定积分int(S,v,a,b),第3讲矩阵、数组和符号运算,a，b为inf时为广义积分,syms x x1 alpha u t;A=cos(x*t),sin(x*t);-sin(x*t
18、),cos(x*t)A=cos(x*t),sin(x*t)-sin(x*t),cos(x*t)int(A,t)ans=1/x*sin(x*t),-cos(x*t)/x cos(x*t)/x,1/x*sin(x*t)int(x1*log(1+x1),0,1)ans=1/4,练习：计算以下积分结果,syms xy=1/(4*x2+2*x+3)z=int(y,-inf,inf)double(z),Ans0.9472,r=symsum(s)%对符号表达式s中的符号变量k（由命令findsym(s)确定的）从0到k-1求和r=symsum(s,v)%对符号表达式s中指定的符号变量v从0到v-1求和r=s
19、ymsum(s,a,b)%对符号表达式s中的符号变量k（由命令findsym(s)确定的）从a到b求和r=symsum(s,v,a,b)%对符号表达式s中指定的符号变量v从a到b求和,syms k n;simplify(symsum(k)ans=1/2*k2-1/2*k simplify(symsum(k2,0,n)ans=1/3*n3+1/2*n2+1/6*n simplify(symsum(k2,0,10)ans=385,c.级数（1）级数求和函数,练习：试求以下级数和,syms kf=k2symsum(f,k,0,10),385,1、积分变换（1）傅立叶变换及其逆变换 B、对函数进行傅立
20、叶逆变换的计算公式为：,7、高级符号运算功能,1、积分变换（1）傅立叶变换及其逆变换A、对函数进行傅立叶变换的计算公式为：,函数：fourier,函数：ifourier,A、傅立叶变换fourier格式：F=fourier(f):对自变量为x的表达式f(x)进行变换，返回值为F()。F=fourier(f,v):指定变换结果为变量v的函数。F=fourier(f,u,v)：指定对函数表达式作关于变量u的傅立叶变换。,F(w)=int(f(x)*exp(-i*w*x),x,-inf,inf),F(v)=int(f(x)*exp(-i*v*x),x,-inf,inf),F(v)=int(f(u)*
21、exp(-i*v*u),u,-inf,inf),例：对以下函数求傅立叶变换,syms xy=exp(-x2)pretty(fourier(y),1/2 2 pi exp(-1/4 w),b、对函数L(s)进行拉普拉斯反变换，计算公式为：,(1)拉普拉斯变换及其反变换A、对函数f(x)进行拉普拉斯变换，计算公式为：,函数：laplace,函数：ilaplace,A、Laplace变换L=laplace(F)x s，F为待进行拉普拉斯变换的代数表达式，默认变量为t。函数返回默认自变量为s。L=laplace(F,z)x z，指定返回结果为自变量z的函数。L=laplace(F,w,z)w z，要求
22、F对自变量w进行拉普拉斯变换，返回值的自变量为z。,第3讲矩阵、数组和符号运算,L(s)=int(F(t)*exp(-s*t),t,0,inf),L(z)=int(F(x)*exp(-z*x),x,0,inf),L(z)=int(F(w)*exp(-z*w),w,0,inf),试对函数 进行拉普拉斯变换，并将所得结果表示为变量v的函数,syms x s w z laplace(sin(x)ans=1/(s2+1)laplace(sin(x),w)ans=1/(w2+1),laplace(sin(x*w)ans=w/(s2+w2)laplace(sin(x*w),w,z)ans=x/(z2+x2
23、),syms x t s vF=sin(x*t+2*t)L1=laplace(F,s,v),L1=sin(x*t+2*t)/v,已知某控制系统的微分方程 其中T0.5s，K10.设初始条件为0，试求系统单位阶跃响应曲线,syms s k tGs=k/(1+t*s)Gs1=subs(Gs,k,t,10,0.5)G=Gs1*1/sy=ilaplace(G)figure,ezplot(y,0 5)xlabel(t)ylabel(y(t)title(系统单位阶跃响应)legend(20*exp(-t)*sinh(t)legend(y(t)=20*exp(-t)*sinh(t),符号运算练习之拉普拉斯应
24、用举例,B、Laplace反变换F=ilaplace(L)L为待进行拉普拉斯反变换的代数表达式，默认变量为s。函数返回默认自变量为t。如果L=L(t),则返回结果为F(x).F=ilaplace(L,y)指定返回F为自变量y的函数。积分对s进行F=ilaplace(L,y,x)要求对自变量y求L的反拉普拉斯变换，返回值的自变量为x。,第3讲矩阵、数组和符号运算,F(t)=int(L(s)*exp(s*t),s,c-i*inf,c+i*inf),F(y)=int(L(y)*exp(s*y),s,c-i*inf,c+i*inf),F(y)=int(L(y)*exp(x*y),y,c-i*inf,c
25、+i*inf),试求以下函数的拉普拉斯反变换，并把L3、L4表示为变量x的函数。,syms s t w x yF1=ilaplace(1/(s-1)F2=ilaplace(1/(t2+1)F3=ilaplace(t(-5/2),x)F4=ilaplace(y/(y2+w2),y,x),F1=exp(t)F2=sin(x)F3=4/3*x(3/2)/pi(1/2)F4=cos(w*x),练习 线性系统冲击响应,编写求任意高阶连续常系数线性系统脉冲响应的程序。解:任意阶次的连续线性系统可用下列微分方程表述：写成传递函数形式为其特性可用系统传递函数的分子分母系数向量b和a来表示。,传递函数反变换的求
26、法,如果分母系数多项式没有重根，则可以把两个多项式之比分解成n个一阶部分分式之和。即：其中p1,p2,pn是分母多项式的n个根，而r1,r2,rn是则对应于这n个根的留数。一阶分式的反变换可以查表得到，容易写出冲击响应的公式如下：可见只要求出根向量p和留数向量r，线性方程的解就得到了。求根是代数问题，当阶次很高时，代数方程没有解析解。可喜的是MATLAB提供了用数值方法求根和留数的函数residue.m，它的调用格式为：r,p=residue(b,a),课堂练习,给出系统的传递函数为求系统的脉冲响应h(t)，如右图所示。,求解1,a=input(多项式分母系数向量a=(取 poly(0,-1,
27、-2,-5)）)b=input(多项式分子系数向量b=(取 1,7,1）)r,p=residue(b,a),%求留数disp(解析式h(t)=r(i)*exp(p(i)*t)k=input(是否要求波形？是,键入1;否,键入0);if k=1dt=input(dt=(取0.05);tf=input(tf=(取5);%设定时间数组 t=0:dt:tf;h=zeros(1,length(t);%h的初始化 for i=1:length(a)-1%根数为a的长度减1 h=h+r(i)*exp(p(i)*t);%叠加各根分量 endFigure,plot(t,h),gridelse,end,求解2,s
28、yms s Gs=(s2+7*s+1)/(s*(s+1)*(s+2)*(s+5)G=Gs*1y=ilaplace(G)figure,ezplot(y,0 5)xlabel(t)ylabel(y(t)title(系统单位脉冲响应),如何求解系统阶跃响应？,求解3,s=tf(s)G=tf(1 7 1,1 8 17 10 0);subplot(3,1,1)impulse(G);%脉冲响应subplot(3,1,2);step(G);%阶跃响应subplot(3,1,3);G1=G/s2;impulse(G1);%斜坡响应,练习：阅读下列程序写出功能,syms s Gs=(s2+7*s+1)/(s*(
29、s+1)*(s+2)*(s+5)G=Gs*1G1=Gs*1/sG2=Gs*1/s2y=ilaplace(G)y1=ilaplace(G1)y2=ilaplace(G2)figure,subplot(3,1,1),ezplot(y,0 7)subplot(3,1,2);ezplot(y1,0 1500)subplot(3,1,3);ezplot(y2,0 300),课堂练习,给出以下程序执行结果S1=sym(1.234)S2=sym(1.234)F3=123fS3=sym(F3)F4=123;456;789S4=sym(F4),给出以下程序执行结果1.234617/500F3=123f Erro
30、r using F4=123;456;789Error in,Sym(f)中，f不能是字符矩阵或者非法的表达式,三、Matlab流程控制结构,分支结构(if)分支结构(switch)for循环while循环嵌套结构程序规范,分支结构(if),if 判断条件1 执行语句1elseif 判断条件2 执行语句2elseif 判断条件3 执行语句3else 执行语句n+1end,只要表达式内所有元素为真，则执行命令串。执行时，依次检查各表达式，只执行第一个表达式为真的命令串，接下来的关系式不检验，跳过其余的if-else-end,分支结构(if),if 判断条件1 执行语句1elseif 判断条件2
31、执行语句2elseif 判断条件3 执行语句3else 执行语句n+1end,if temp 10,disp(cold!);elseif temp 20,disp(cool);elseif temp 27,disp(warm);elseif temp 31,disp(hot);else disp(too hot);end,分支结构(switch),switch 表达式case 表达式1语句段1case 表达式2,表达式3,表达式m语句段2otherwise语句段nend,sco=input(Please input the scoren);switch floor(sco/10)case 9,
32、10 t=A;case 8 t=B;case 7 t=C;case 6 t=D;case 0,1,2,3,4,5 t=E;otherwise disp(Invalid score)end disp(the core is t),输入一个百分制成绩，要求输出成绩等级A、B、C、D、E。其中：90100分为A,8089分为B，7079分为C，6069分为D，60分一下为E。要求：用switch语句实现，,for循环,循环结构for end循环格式1:for 循环变量=初值:增量:终值语句段end 格式2:for 循环变量=向量语句段end,for循环举例1,for 循环变量=循环值 循环操作end
33、,s=0;for k=1:100,s=s+1/k2;end,for k=hello world,fprintf(%c,k);pause(0.1)end,for k=magic(3),k+1end,向量功能,for循环练习,两种格式：求S=1+2+100 S=0;for x=1:1:100S=S+x;end S S=0;v=1:100;for x=vS=S+x;end S,while循环,while 循环条件（逻辑表达式）循环操作end while 1,if 循环条件,break;end 循环操作end,n=input(n=);while n=1,if mod(n,2)=1,n=n*3+1;br
34、eak;else n=n/2;break;endend,表达式内的元素为真,练习：,用循环结构求解,s=0;for i=i+100s=s+i;end,s=0;i=1;while i=100s=s+i;i=i+1;end,Sum(1:100),用满足 的最小m值,s=0;m=0;while(s=10000)m=m+1;s=s+m;end,求解级数求和,tici=1:100000;s=sum(1./2.i+1./3.i);toc,tics=0;for i=100000 s=s+1/2i+1/3iendtoc,比较速度,y=0;i=1;while 1f=1/(2*i-1);y=y+f;if y3br
35、eak;endi=i+1end,求y3时的最大n的值，并求此时n值和n值对应的y值。,M文件M文件有两种格式命令文件（script file），另外这类程序包含一连串的命令，执行时依次进行。函数文件（function file）第一句可执行语句是以function引导的定义语句，函数文件中的变量都是局部变量，仅在函数内部有效。函数名和文件名必须相同,程序调试：Debug breakpoint,Step：StepinStep outRunContinueGo until cursorExit debug mode,set/clear breakpointClear all breakpointS
36、top if errorStop if warningStop if all error,程序规范,程序按照结构缩进(matlab会自行处理)按照逻辑关系分段添加适当的注释(5%-30%)规范所有的标识符(文件名、变量名)适当增加空行、空格以保证可读性循环语句在matlab中执行速度较慢，所以能用矩阵运算时，尽量不要用循环结构，提高代码效率办法:参考Matlab内部程序风格,编程要点,尽量避免使用循环，而使用向量或矩阵如果要使用循环，在循环语句前也要尽量对向量、矩阵或数组预先用ones或zeros进行内存分配尽量使用matlab内部函数或工具箱函数。绝大多数常见数学计算都可在matlab中找到
37、相应的函数命令,（1）设，对f(x)和g(x)进行因式分解。（2）运用 simplify 命令简化，并与采用 simple 命令得出的结果相比较。（3）计算下列极限 a.b.（4）建立一个符号表达式 y=sin(a+b)2*xc)，分别以变量 c 进行积分、以变量 x 从pi/2 到 pi 进行积分，（5）求符号矩阵 F 对变量 x 的一阶微分、对变量 a 的二阶微分 F=cos(a*x)sin(a+b)2*x)-sin(b*x)en;,第3讲上机练习（5）,（6）建立一个包含 的符号表达式，并将该表达式转换为符号变量。（7）用符号计算验证下述三角恒等式的正确性。（8）用 rand 函数建立一
38、个 33 维随机数值矩阵，并将该数值矩阵转换为符号矩阵，比较它们的不同。（9）分别用 1 位和 10 位有效数字计算 2*pi+3.1415 的值。（10）求下述矩阵的逆、特征根和矩阵的行列式值。,上机练习（5）,（11）编写程序：在2-100中，找出不能被3，7，11整除的数，输出第2个数（12）阅读下列程序，并将其编辑成m文件，命名；用debug下的断点（breakpoint）设定、取消，单步运行的工具进行调试，并画出流程图 i=1;for n=2:100 if mod(n,2)=0|mod(n,3)=0|mod(n,5)=0 continue else X(i)=n;i=i+1;endend,
