矩阵、数组和符号运算.ppt

第3讲 矩阵、数组和符号运算2 第3讲 Matlab流程控制结构3 第3讲矩阵、数组和符号运算 掌握内容： （1）了解 MATLAB 的符号变量，掌握 MATLAB 符号表 达式、符号矩阵的两种创建方法。 （2）掌握 MATLAB 符号数学函数的创建。 （3）掌握符号矩阵的基本运算及MATLAB 关于不同精度 的控制方法。 （4）掌握符号微积分内容，包括求函数的极限、对符号表 达式求导数和微分、符号积分、符号求和、傅立叶变换及 其逆变换等。 （5）掌握各种符号方程的求解方法和函数命令。 （6）了解 MATLAB 可视化的符号函数分析界面及使用。 （7）初步了解 MAPLE 的符号资源。 二、符号及运算 第3讲矩阵、数组和符号运算 抽象运算：公式推导、因式分解、求解代数方程或微分方 程的精确解 符号数学工具箱 1）通过基本符号数学工具箱的专用函数； 符号表达式和符号矩阵的操作； 多项式的化简、展开和代入； 线性代数； 微积分； 符号方程的求解； 特殊的数学函数。 2）通过 maple.m、mpa.m 两个专门设计的 M 文件进行符 号运算； 3) 通过 MATLAB 中的函数计算器（Function Caculator ）。 第3讲矩阵、数组和符号运算 1、符号变量的创建 a.符号矩阵定义函数 sym 函数 S=sym(arg) ，从表达式 arg 创建一个 sym 对 象 S x=sym(x) x = sym(x,real) x = sym(x,unreal) 附加属性 pi = sym(pi) delta = sym(1/10) 第3讲矩阵、数组和符号运算 A=2/5,4/0.78,sqrt(23)/3;0.33,0.3333,log(4) 输入数值矩阵A A = 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) 33/100, 3333/10000, 6243314768165359*2(-52) 不管数值矩阵的元素是以分数或是浮点数表示，转换后 的符号矩阵都将以最接近有理式的形式给出。 b.符号定义函数 syms 函数 syms arg1 arg2 . syms a b c x y 2、符号表达式和矩阵的创建 a.字符串直接输入创建 符号表达式和符号方程对空格很敏感。因此，在创建符 号表达式或符号方程时，不要在字符间任意加空格符； 符号计算中出现的数字也是当作符号处理的； f=a*x2+b*x+c f = a*x2+b*x+c f=a*x2+b*x+c=0 f = a*x2+b*x+c=0 第3讲矩阵、数组和符号运算 第3讲矩阵、数组和符号运算 这种方法输入符号矩阵与字符串矩阵的输入相似。但要 保证在同一列中各元素字符串有同样的长度，在较短的字 符串前后用空格符填充； 这种方法要求符号矩阵每一行的两端都有方括号，而字 符串矩阵仅在首尾有方括号。 B=4+x x2 x ;x3 5*x-3 x*a B = 4+x x2 x x3 5*x-3 x*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+c f =a*x2+b*x+c syms x a B=4+x x2 x;x3 5*x-3 x*a B = 4+x, x2, x x3, 5*x-3, x*a 不能创建符号方程 符号方程：equ1=(a*x+b=0) a. 基本运算 四则运算 两个符号矩阵的大小相等方可进行加减运算，符号 矩阵和符号标量的加减运算按照数组运算规则进行 ； 两个符号矩阵只有内积相等时才可以进行乘法运算 ； 符号的乘方运算 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-a ans = x-1/x, 1-1/(x+1) x+2-1/(x+2), -1/(x+3) ab ans = -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.b ans = x2, x+1 (x+2)2, 0 a2 ans = 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） 因式分解法（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+(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*t E2 = cos(x)*cos(y)-sin(x)*sin(y) E3 = 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 y H1 = 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)*x H2 = x2+x*y； (-2+y2)*y （4）因式分解法（factor） 格式 factor(s) %符号表达式s的因 式分解函数，说明： S为符号矩阵或符 号表达式，常用于多项式的因式分解。 将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 c R1 = 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)2 R2 = (a+b)(1/2*c) 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（sincos）：将代数式尽量转化为由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 y n,d=numden(x/y+y/x) n = x2+y2 d = 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) n1 =6461369247334093 d1 =9007199254740992 n2 = x2+y2 d2 = y*x n3 = a, 1； 1, d d3 = 1, b； c, 1 C.符号变量（表达式）的提取与代入 （1）subexpr 格式 Y,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 y subs(x+y,x,4) ans =4+y syms a b f=subs(cos(a)+sin(b),a,b,sym(alph a),2) f =cos(alpha)+sin(2) 试将以下表达式中的变量x分别用y和m nt代替。 f=a*x2+b*x+c syms x y a b c f=a*x2+b*x+c sf1=subs(f,x,y) sf2=subs(f,x,m+n*t) a. 符号极限 limit(F,x,a)：求当xa时表达式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) limit(1/x,x,0,right) ans = inf 6、符号微积分 Matlab自变量确定原则：除i、j外，字母位置最接近x的小 写字母为自变量；如果表达式中没有变量，x会被视为默认 的变量。由函数findsym可一找到默认变量 练习：求下列函数极限 syms x h f=(sin(x+h)-sin(x)/h limit(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),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 x y=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 = symsum(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 k f=k2 symsum(f,k,0,10) 385 1、积分变换 （1）傅立叶变换及其逆变 换 B、对函数进行傅立叶逆 变换的计算公式为： 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)*exp(-i*v*u),u,-inf,inf) 例：对以下函数求傅立叶变换 syms x y=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，要求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) syms x t s v F=sin(x*t+2*t) L1=laplace(F,s,v) L1 = sin(x*t+2*t)/v 已知某控制系统的微分方程 其中T0.5s，K10.设初始条件为0，试求系统单 位阶跃响应曲线 syms s k t Gs=k/(1+t*s) Gs1=subs(Gs,k,t,10,0.5) G=Gs1*1/s y=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) 符号运算练习之拉普拉斯应用举例 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+i*inf) 试求以下函数的拉普拉斯反变换，并把L3、L4 表示为变量x的函数。 syms s t w x y F1=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来表 示。 传递函数反变换的求法 如果分母系数多项式没有重根，则可以把两个多项式 之比分解成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,-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=1 dt=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); % 叠加各根分量 end Figure,plot(t,h),grid else, end 求解2 syms s Gs=(s2+7*s+1)/(s*(s+1)*(s+2)*(s+5) ) G=Gs*1 y=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*(s+1)*(s+2)*(s+5) G=Gs*1 G1=Gs*1/s G2=Gs*1/s2 y=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=123f S3=sym(F3) F4=123;456 ;789 S4=sym(F4) 给出以下程序执行 结果 1.234 617/500 F3=123f Error using F4=123;456; 789 Error in Sym(f)中，f不能是字符矩阵或者非法的表达式 三、Matlab流程控制结构 分支结构(if) 分支结构(switch) for循环 while循环 嵌套结构 程序规范 分支结构(if) if 判断条件1 执行语句1 elseif 判断条件2 执行语句2 elseif 判断条件3 执行语句3 else 执行语句n+1 end 只要表达式内所有元素为 真，则执行命令串。 执行时，依次检查各表达 式，只执行第一个表 达式为真的命令串， 接下来的关系式不检 验，跳过其余的if- else-end 分支结构(if) if 判断条件1 执行语句1 elseif 判断条件2 执行语句2 elseif 判断条件3 执行语句3 else 执行语句n+1 end if temp S = 0; for x=1:1:100 S = S + x; end S S = 0; v = 1:100; for x = v S = 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; break; else n = n / 2; break; end end 表达式内的元素为真 练习： 用循环结构求解 s=0; for i=i+100 s=s+i; end s=0;i=1; while i3 break; end i=i+1 end 求求y3y3时的最大时的最大n n的的 值，并求此时值，并求此时n n值值 和和n n值对应的值对应的y y值值 。 M文件 M文件有两种格式 命令文件（script file），另外这类程 序包含一连串的命令，执行时依次进行 。 函数文件（function file）第一句可 执行语句是以function引导的定义语 句，函数文件中的变量都是局部变量， 仅在函数内部有效。函数名和文件名必 须相同 程序调试： Debug breakpoint StepStep： StepinStepin Step outStep out Run