1、1第三讲 MATLAB 的符号运算（注：文中红色字体为命令执行的结果，在 Command 窗口中显示）3-1 符号对象的创建和使用1符号运算入门 符号运算的特点是，运算过程中允许存在非数值的符号变量。先看如下示例：函数 ，用 MATLAB 求它的微积分，命令如下：2)(sin)xff=sin(x)2; %定义符号函数 f(x)dfdx=diff(f) %求 的指令dxf)(intf=int(f) %求 的指令显示的计算结果为：dfdx=2*sin(x)*cos(x)intf=-1/2sin(x)*cos(x)+1/2*x所以， ， 。此例中，首先定义符号函xdxfcosin)(2 xxdfco
2、sin)(21数 f=sin(x)2，然后由符号运算获得 的微分和积分。)si2定义符号变量 在使用符号变量之前，应先声明某些要用到的变量是“符号”变量。声明符号变量的语句：syms 变量名列表或： sym(变量名)其中各个变量名应该用空格分隔，而不能用逗号分隔。如创建符号变量 x 和 a：x=sym(x)a=sym(alpha)或用： syms x a %定义符号变量 x 和 a这里，变量 x 和 a 的类型是符号对象，它们被定义后，即可参与符号运算。3定义符号表达式和符号方程 符号表达式和符号方程是两种不同的操作对象。区别在于：符号表达式不包含等号（=） ，而符号方程须带等号。它们的创建方
3、式相同。如：要考虑二次函数 f=ax2+bx+c，可以创建符号表达式，赋值给符号变量 f。f=sym(a*x2+b*x+c)2或： f=a*x2+b*x+c此例中，将符号表达式赋给符号变量 f，但这不是必需的，引入符号变量是为了以后调用方便。在这种情况下，没有创建对应于表达式中 a、b、c、x 项的变量，为了执行符号数学运算（如微分、积分等） ，必须显式地创建这些变量，可用下列命令创建：syms a b c x如下例中创建了符号表达式和符号方程，分别赋给相应的符号对象。syms x a b cf=sin(x)2; %创建符号表达式 sin(x)2 赋给变量 feq=a*x2+b*x+c=0 %
4、创建的符号方程赋给变量 eq4定义抽象函数和符号数学函数若要创建抽象函数 f(x)，可使用：f=sym(f(x)则 f 就象 f(x)一样参与运算。在 Symbolic Math Toolbox 中可利用符号表达式创建符号数学函数。如下面的命令将产生符号表达式 r、t、f ：syms x y zr=sqrt(x2+y2+z2) %或 r=sym(sqrt(x2+y2+z2)t=atan(y/x)f=sin(x*y)/(x*y)3-2 数值与符号的转换Sym 函数有四种选项将数值结果转换为符号表达式。如对于变量 rho，定义为： ，在 MATLAB 命令窗口中定义 rho:251rhorho=(
5、1+sqrt(5)/2将这一数值结果转化为符号表达式：（1）sym(rho,f) 返回符号浮点表示形式，结果为：sym(rho,f)ans=1.9e3779b97f4a8*2(0)（2）sym(rho,r) 返回符号有理数表示形式，这是 sym 的默认设置，当调用 sym 而没有第二个参数时，相当于 sym 使用了 r 选项。结果为：sym(rho,r) ans=7286977268806824*2(-52)3例如：sym(0.75)ans=3/4（3）sym(rho,e) 返回符号有理数表示形式，同时根据 eps 给出 rho 的理论表达式和实际计算的差。sym(rho,e)ans=7286
6、977268806824*2(-52)这里，rho 的理论值和实际浮点值相同，对于 1/3，可得：sym(1/3,e)ans=1/3-eps/12（4）sym(rho,d) 返回符号十进制小数表示形式。有效位数由 digits 定义。Digits 的默认值是 32 位。sym(rho,d)ans=1.6180339887498949025257388711907如果想要一个较短的结果，可由 digits 定义有效位数。digits(5)sym(rho,d)ans=1.6180上述四种格式中，最后得到的都是符号对象变量。3-3 符号算术运算1定义符号矩阵MATLAB 中的数值矩阵不能直接参与符号
7、运算，必须经过转换。 Sym 函数的一个非常有用的功能就是将数值矩阵转换成符号矩阵。例如，下面的代码将产生一个 33 的 Hilbert 矩阵。A=hilb(3)A=1.0000 0.5000 0.33330.5000 0.3333 0.25000.3333 0.2500 0.2000对 A 使用 sym 命令可获得一个具有无穷精度的 33 的符号形式的矩。A=sym(A)A= 1,1/2,1/31/2,1/3,1/41/3,1/4,1/54当调用 Symbolic Math Toolbox 时，默认地调用有理算术运算。当使用 sym 函数创建符号变量时，有理算术运算将被启动。例如，对于双精度
8、矩阵 A，sym 函数将转化为对应的符号矩阵形式。A=1.1000 1.2000 1.30002.1000 2.2000 2.30003.1000 3.2000 3.3000A=sym(A)A=11/10, 6/5,13/1021/10,11/5,23/1031/10,16/5,33/10对于非数值的符号矩阵，通过 sym 函数，可用与普通矩阵定义类似的格式来定义。如，下例将定义符号矩阵 A、符号矩阵 B、符号矩阵 C：Sym a b c d xA=sym(a b;c d)A=a,bc,dB=sym(2*a+b,3*b;5*c+d,2*d)B=2*a+b, 3*b5*c+d, 2*dC=sym
9、(1/(1+a),sin(x),(b-x)/(a+x);1,exp(x),x2)C=1/(1+a), sin(x), (b-x)/(a+x) 1, exp(x), x22符号矩阵的加、减运算实现符号矩阵加、减的指令是：symadd(A,B) %给出两个符号矩阵的和 A+Bsymsub(A,B) %给出两个符号矩阵的差 A-B显然，参与符号加、减运算的符号矩阵应具有相同的维数。如，对于上面定义的符号矩阵 A 和 B，利用加、减指令的结果是：symadd(A,B)ans=3*a+b, 4*b6*c+d, 3*dsymsub(A,B)ans=5 -a-b, -2*b-4*c-d, -d对于具有数值的
10、符号矩阵的加法，使用同样计算方法。如计算两个符号矩阵的和：A=hilb(3)A=1.0000 0.5000 0.33330.5000 0.3333 0.25000.3333 0.2500 0.2000A=sym(A)A= 1, 1/2, 1/3 1/2, 1/3, 1/4 1/3, 1/4, 1/5B=sym(2*a,b,c;2*b,c;a,b,2*c)B=2*a, b, c a, 2*b, c a, b, 2*csymadd(A,B)ans=1+2*a, 1/2+b, 1/3+c 1/2+a, 1/3+2*b, 1/4+c 1/3+a, 1/4+b, 1/5+2*c3符号矩阵的乘、除运算实现
11、符号矩阵乘、除的指令是：symmul(A,B) %给出两个符号矩阵的积 A*Bsymdiv(A,B) %给出两个符号矩阵的商 A/B如，对于上面定义的符号矩阵 A 和 B，利用乘、除指令的结果是：symmul (A,B)ans=a*(2*a+b)+b*(5*c+d), 3*a*b+2*b*dc*(2*a+b)+d*(5*c+d), 3*c*b+2*b2symdiv(A,B)ans=(2*a*d-b*d-5*c*d)/(-15*c*b-b*d+4*a*d), -b*(-b+a)/(-15*c*b-b*d+4*a*d) -d*(3*c+d)/ (-15*c*b-b*d+4*a*d), (2*a*d
12、-3*c*b-b*d) /(-15*c*b-b*d+4*a*d)4符号变量替换有时，一个符号解或一个符号表达式中，需将一些符号变量替换成数字或其它符号，可利用 subs 函数实现。6subs 函数适用于单个符号矩阵、符号表达式、符号代数方程和微分方程。subs 函数有两种格式：（1）subs(S,NEW)用新变量 NEW 替换 S 中的默认变量。如：符号矩阵 G，用pi/3替换G 中的默认变量 x：G=sym(a*sin(b+x),a+b,exp(a*x),sqrt(x);G1=subs(G,pi/3)G1=a*sin(b+ pi/3), a+b, exp(a* pi/3), sqrt(pi/
13、3)和利用 pi/3 替换 G 中的默认变量 x 比较：G=sym(a*sin(b+x),a+b,exp(a*x),sqrt(x);G1=subs(G,pi/3)G1=a*sin(b+1/3* pi), a+b, exp(1/3*a* pi), 1/3*3(1/2)*pi(1/2)（2）subs(S,NEW,OLD)用新变量 NEW 替换 S 中的指定的变量 OLD。如：符号表达式f，用 2 替换 f 中的变量 a：f=sym(sin(1/3*a*pi);subs(f,2,a)ans=sin(1/3*2*pi)如果用数值替换符号表达式中的变量，将得到相应变量由数值替代的表达式或数值。如：f=s
14、ym(sin(1/3*a*pi);subs(f,2,a)ans=0.86603-4 符号微积分运算1确定符号变量当进行数学运算时，对应变量的选取很容易得到。如，对于函数 f=xn，当对 f 求导时，自然地是对 x 求导，n 看成常数。而在 MATLAB 中，如何知道是对 x 求导而不是对 n 求导呢？它通过符号表达式中隐含的符号变量来确定。在 Symbolic Math Toolboxs 中，确定一个符号表达式中的符号变量的规则是：（1） 只对（除 i,j 外）单个小写英文字母进行检索。（2） 小写字母 x 是首选符号变量。（3） 其余小写字母被选用符号变量的次序：在英文字母表中，靠近“x”的
15、优先，在“x”之后的优先。按照这一规则，对 f=xn 求导时，自然是对 x 求导，n 看成常数。7工具箱中还提供了 findsym 函数来确定表达式中的符号变量。如， findsym(f,1)寻找第一个符号变量。findsym 函数中的第二个参数，代表了在符号对象中想要寻找的符号变量的个数。默认时，将给出符号表达式中的所有符号变量。如下例给出了 f 中的符号变量。syms a b c xf=sym(a*x2+b*x+c);findsym(f,1)ans=xfindsym(f,2)ans=x,c2符号微分运算对符号表达式微分的函数是 diff()。该函数可以求符号表达式的一阶导数、n 阶导数。调
16、用格式：（1）diff(f) %传回 f 对变量 x 的一次微分值 syms a x;f=sin(a*x);df=diff(f)df=cos(a*x)*a或：f=sin(a*x)df=diff(f)df=cos(a*x)*a（2）diff(f,a) %传回 f 对指定变量 a 的一次微分值 下列命令分别计算 f 对变量 x 和 n 的微分：syms x nf=xndiff(f,x)ans=xn*n/xdiff(f,n)ans=xn*log(x)（3）diff(f,n) %传回 f 变量 x 的 n 次微分值 diff(f,a,n) %传回 f 对指定变量 a 的 n 次微分值 syms a x
17、;f=sin(a*x);df=diff(f,x,2)8df=-sin(a*x)*a2diff 函数也可以使用符号矩阵作为它的输入，此时，微分按矩阵元素逐个进行。数值微分函数也是用 diff，因此这个函数是靠输入的参数决定是以数值或是符号微分，如果参数为向量则执行数值微分，如果参数为符号表示式则执行符号微分。 先定义下列三个方程式，接着再演算其微分项： S1 = 6*x3-4*x2+b*x-5; S2 = sin(a); S3 = (1 - t3)/(1 + t4); diff(S1) ans=18*x2-8*x+b diff(S1,2) ans= 36*x-8 diff(S1,b) ans=x
18、 diff(S2) ans=cos(a) diff(S3) ans=-3*t2/(1+t4)-4*(1-t3)/(1+t4)2*t3 simplify(diff(S3) ans=t2*(-3+t4-4*t)/(1+t4)2 3符号积分运算int 函数用以计算函数的积分项，这个函数要找出一符号式 F 使得 diff(F)=f。如果积 分式的解析式 (analytical form, closed form) 不存在的话或是 MATLAB 无法找到，则 int 传回原输入的符号式。相关的函数语法有：（1）int(f) %传回符号表达式 f 对变量 x 的不定积分 syms x;s=1/(1+x2)
19、;f=int(f)f=atan(x)（2）int(f,t) %传回符号表达式 f 对指定变量 t 的不定积分syms a;9s=sin(a*u);int(s,a)ans=-1/u*cos(a*u)（3）int(f,a,b) %传回符号表达式 f 对变量 x 的积分值，积分区间为a,b，a 和 b 为数值式 int(f,t,a,b) %传回符号表达式 f 对变量 t 的积分值，积分区间为 a,b，a 和 b 为数值式 int(f,m,n) %传回 f 对变量 x 的积分值，积分区间为 m,n，m 和 n 为符号式 我们示范几个例子： S1 = 6*x3-4*x2+b*x-5; S2 = sin(
20、a); S3 = sqrt(x); int(S1) ans=3/2*x4-4/3*x3+1/2*b*x2-5*x int(S2) ans= -cos(a) int(S3) ans= 2/3*x(3/2) int(S3,a,b) ans=2/3*b(3/2)- 2/3*a(3/2) int(S3,0.5,0.6) ans=2/25*15(1/2)-1/6*2(1/2) 4符号微积分运算示例例 1：对于函数 ，求 。yxfsin),(yxf),(syms x,y;f=x*siny;dfdxdy=diff(diff(f,x),y)dfdxdy=cos(y)例 2：对于函数 ，先求 s 关于 x 的不
21、定积分，再求所得结果关于 y 的不定积分，xyes),(即计算 。dxey(syms x,y;s=x*exp(-x*y);int(int(s,x),y);10ans=1/y*exp(-x*y)即 。yedxexy)(例 3：给定一个函数： ，求该函数的微分，然后再将微分结果进行积分运算。sin3412syms x y y1 y2y=1/(x2+4*x+3)*sinx;y1=diff(y,x)y1=-1/(x2+4*x+3)2*sinx*(2*x+4)y2=int(y1,x)y2=1/(x2+4*x+3)*sinx3-5 符号绘图函数 ezplot()例：创建符号函数 ，并绘该函数的波形。xfc
22、os)(451x=sym(x);f=sym(1/(5+4*cos(x);ezplot(f)-6 -4 -2 0 2 4 60.10.20.30.40.50.60.70.80.91x1/(5+4 cos(x)求 f 的微分并绘出结果的波形f1=diff(f)f1 =4/(5+4*cos(x)2*sin(x)ezplot(f1)11-6 -4 -2 0 2 4 6-1-0.8-0.6-0.4-0.200.20.40.60.81x4/(5+4 cos(x)2 sin(x)求 f 的二次微分并绘出结果的波形f2=diff(f,2)f2 =32/(5+4*cos(x)3*sin(x)2+4/(5+4*c
23、os(x)2*cos(x)ezplot(f2)-6 -4 -2 0 2 4 6-1-0.8-0.6-0.4-0.200.20.40.60.81x32/(5+4 cos(x)3 sin(x)2+4/(5+4 cos(x)2 cos(x)求 f2 的二次积分并绘出结果的波形g=int(int(f2)g =-8/(tan(1/2*x)2+9)ezplot(g)-6 -4 -2 0 2 4 6-0.9-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.10x-8/(tan(1/2 x)2+9)此时，我们发现，g 和原函数 f 形式不同，但图形相似。下面求 g 和 f 的差值：r=f-gr =
24、1/(5+4*cos(x)+8/(tan(1/2*x)2+9)12ezplot(r)-6 -4 -2 0 2 4 611111x1/(5+4 cos(x)+8/(tan(1/2 x)2+9)可见，图形窗口中绘制的 f-g 的波形是一条幅度为 1 的直线。对 r 进行简化，e=simple(r)e=1所以，f-g 为是个恒定常数。3-6 求解常微分方程式 MATLAB 解常微分方程式的语法是 dsolve(equation,condition)，其中 equation 代表常微分方程式即 y=g(x,y)，且须以 Dy 代表一阶微分项 y D2y 代表二阶微分项 y，condition 则为初始
25、条件。 假设有以下三个一阶常微分方程式和其初始条件 y=3x2，y(2)=0.5 y=2xcos(y)2，y(0)=0.25 y=3y+exp(2x)，y(0)=3 对应上述常微分方程式的符号运算式为： soln_1 = dsolve(Dy = 3*x2,y(2)=0.5) ans= x3-7.500000000000000 ezplot(soln_1,2,4) % 看看这个函数的长相 soln_2 = dsolve(Dy = 2*x*cos(y)2,y(0) = pi/4) ans= atan(x2+1) soln_3 = dsolve(Dy = 3*y + exp(2*x), y(0) =
26、 3) ans= -exp(2*x)+4*exp(3*x) 133-7 连续信号的微分与积分连续信号的微分可由 diff 近似计算：h=0.001;x=0:h:pi;y=diff(sin(x.2)/h;定积分可由 quad 函数或 quadl 函数实现:quad(function_name,a,b) % function_name 为被积函数名，a 和 b 指积分区间例：求 example2-5(第二讲中) 所示三角波求微分和积分第二讲中三角波的产生：%example2-5 t=-4:0.001:4;yt=tripuls (t,4,0.5);plot(t,yt) 先将三角波写成 MATLAB
27、函数（m 文件） ，名为 x2_5:function yt=x2_5(t)yt=tripuls (t,4,0.5);用 diff 和 quad 函数可实现 三角波的微分和积分%example3-1 differentiation（微分）h=0.001;x=-3:h:3;y1=diff(x2_5(t)*1/h;plot(t(1:length(t)-1),y1)title(dx(t)/dt)%example3-2 integration（积分）t=-3:0.1:3;for n=1:length(t)y2(n)=quad(x1_2 ,-3,t(n);endplot(t,y2)title(integr
28、al of x(t) )143-8 练习1、用符号表达式表示信号 ，并绘制波形。)sin()ttf42、实现下列连续时间信号 x(t)的微分和积分。20135t)(tx3、信号 ，用 MATLAB 符号运算的相关命令求： 、)()()22tuttf )(2tf、 、 、 ，并绘出其时域波形。(2tttftf3 答案：syms t %定义符号变量 tf=sym(t/2+1)*(heaviside(t+2)-heaviside(t-2) %创建符号表达式 f，heaviside 函数见第二讲subplot(2,3,1), ezplot(f,-3,3) %subplot 把图形窗口分成 2 行 3
29、列的 6 个小窗口，当前波形画在第 1 个小窗口中，下同y1=subs(f,t,t+2) %subs()函数指将连续信号 f 中的时间变量 t 用 t+2 替换，下同subplot(2,3,2), ezplot(y1,-5,1)y2=subs(f,t,t-2)subplot(2,3,3), ezplot(y2,-1,5)y3=subs(f,t,-t)subplot(2,3,4), ezplot(y3,-3,3)y4=subs(f,t,2*t)subplot(2,3,5), ezplot(y4,-2,2)y5=-fsubplot(2,3,6), ezplot(y5,-3,3)请注意：上述命令执行后，f、y1、y2、y3、y4、y5 的符号表达式（在 Command 窗口中显15示） 。
