1、1第 2 章 符号计算符号计算：解算数学表达式、方程不是在离散化的数值点上进行，而是凭借一系列恒等式，数学定理，通过推理和演绎，获得解析结果。符号计算建立在数值完全准确表达和推演严格解析的基础之上，所得结果完全准确。特点：一，相对于 MATLAB 的数值计算 “引擎”和“函数库”而言，符号计算的“引擎”和“函数库”是独立的。二，在相当一些场合，符号计算解算问题的指令和过程，显得比数值计算更自然、更简明。三，大多数理工科的本科学生在学过高等数学和其他专业基础课以后，比较习惯符号计算的解题理念和模式。运算引擎 MuPADMuPAD 是极佳数学及符号数值运算绘图软件，同时也作为 MATLAB7.8
2、的符号计算工具箱，是一具有人工智能的数学软件,非常适合科学家及工程师使 用.更适合每一个人使用,使用的方法非常简单,只要输入方程式就立刻得到答案,可以求 Symbolic 符号解,多项式之根,求非线性方程式之根,矩阵及向量 Vector and Matrices 运算,代数 Algebra 运算,求积分之值,求微分之值 Calculus 微积分等。方程式可以处理复数计算 . 完美的绘图功能,图型输入,输出,轻松无比的绘图,可以输入多个 2-D 函数或极坐标函数或 3-D 函数 ,选择所要绘图参数,就可以完成图形，以及图形的动画 制作也是非常方便。数值计算结果并不是 MATLAB 命令行窗口所得
3、的类似代码形式，而是规范数学格式。并拥有一内建的程序语言，帮助文档以及文本操作，文 本操作在一定程度上可以取代 word. 是一个超级的工程数学计算器 .MathWorks 自从 2008 年 10 开始，在 Matlab 的新版本(Matlab2008a，即 7.6 之后)中使用MuPAD 内核替换原来的 Maple 符号计算内核！2.1 符号对象和符号表达式MATLAB 依靠基本符号对象（包括数字、参数、变量 ）、运算符及一些预定义函数来构造和衍生符号表达式和符号方程。2.1.1 符号对象的创建和衍生一 生成符号对象的基本规则 任何基本符号对象都必须借助专门的符号函数指令 sym 或 sy
4、ms 定义。 任何包含符号对象的表达式或方程，将继承符号对象的属性。2二 符号数字符号（类）数字的定义：sym(Num) 创建一个符号数字 Numsc=sym(Num) 创建一个符号常数 sc，该常数值准确等于 Num说明：Num 代表一个具体的数字Num 必须处于（英文状态下的）单引号内，构成字符串（关于字符串参见附录A.1）。【例 2.1-1】符号（类）数字与数值（类）数字之间的 差异。a=pi+sqrt(5) % 创建方式sa=sym(pi+sqrt(5) Ca=class(a) % 类别判断Csa=class(sa) vpa(sa-a) digits digits(8) 三 符号参数表
5、达式 e-axsinbx 中的 a,b 称为参数。定义格式：syms Para 定义符号参数 ParaPara=sym( Para) syms Para Flag 定义具有 Flag 指定属性的符号参数 ParaPara=sym( Para, Flag) syms Para1 Para2 ParaN 定义 Para1 Para2 ParaN 为符号参数syms Para1 Para2 ParaN Flag 定义 Para1 Para2 ParaN 为具有 Flag 指定属性的符号参数 符号参数名不要用处于“字母表中小写字母 x 及其两侧的英文字母”开头。 Flag 表示参数属性，可具体取以下词
6、条：positive 表示那些符号参数取正实数；real 表示那些符号参数限定为实时；unreal /clear 表示那些符号参数为不限定的复数。四 符号变量e-axsinbx 中的 x 称为变量，符号变量的定义同符号参数。确定自由符号变量的规则：3 在专门指定变量名的符号运算中，解题一定围绕指定变量名进行。 自动识别符号变量时，字母的优先次序为 x,y,w,z,v 等。自动识别表达式中自由、独立的符号变量的指令： findsym(EXPR) 确认表达式 EXPR 中所有自由符号变量findsym(EXPR, N) 确认表达式 EXPR 中距离 x 最近的 N 个自由符号变量【例 2.1-2】
7、用符号计算研究方程 的解。02wvzu（1）不指定变量情况syms u v w z % 定义符号参数/变量Eq=u*z2+v*z+w;result_1=solve(Eq) %findsym(Eq,1) （2）指定变量情况result_2=solve(Eq,z) 【例 2.1-3】对独立自由符号变量的自动辨认。（1）syms a b x X Y % 定义符号参数/ 变量k=sym(3); % 符号常数z=sym(c*sqrt(delta)+y*sin(theta1); % 直接定义符号表达式EXPR=a*z*X+(b*x2+k)*Y; % 构成衍生符号表达式（2）findsym(EXPR)（3）
8、findsym(EXPR,1) （4）findsym(EXPR,2),findsym(EXPR,9) 【例 2.1-4】findsym 确定自由变量是对整个矩阵进行的。syms a b t u v x yA=a+b*x,sin(t)+u;x*exp(-t),log(y)+vfindsym(A,5) A = a + b*x, u + sin(t) x/exp(t), v + log(y)ans =x,y,v,u,t 2.1.2 符号计算中的算符 与数值计算中的算符在形状、名称和使用方法上几乎完全相同。 仅注意：在符号对象的关系运算符中，只有算符“=”，“=”4比较结果为“真”时，用 1 表示;
9、否则用 0 表示。2.1.3 符号计算中的函数指令表 2.1-1 MATLAB 中可调用的符号计算函数指令类 别 情况描述 与数值计算对应关系三角函数、双曲函数及反函数；除 atan2 外 名称和使用方法相同指数、对数函数（如 exp，expm） 名称和使用方法相同复数函数（注意：没有幅角函数 angle） 名称和使用方法相同矩阵分解函数（如 eig，svd） 名称和使用方法相同方程求解函数 solve 不同微积分函数(如 diff，int) 不完全相同积分变换和反变换函数（如 laplace，ilaplace） 只有离散 Fourier 变换基本函数绘图函数（如 ezplot，ezsurf）
10、 数值绘图指令更丰富经典特殊函数 如误差函数 erf、贝塞尔函数 besselj、第一类完全椭圆积分 EllipticK 等；通过 mfunlist 可以看到所有经典函数名mfun functions perform numerical, not symbolic, calculations. The input parameters should be scalars, vectors, or matrices of type double, or complex doubles, not symbolic variables.部分Maple 库函数 Maple 库函数在符号计算的扩展目录上
11、；可通过mhelp index 看到各子函数库的名称；函数的数量很大；使用库函数，需要具备 Maple 语言知识注意：使用函数注意数据类型。就数字而言，有双精度和符号类数字之分。2.1.4 符号对象的识别为了函数指令与数据对象的适配，MATLAB 提供了用于识别数据对象属性的指令：class(var) 给出变量 var 的数据类别（如 double，sym 等）isa(var, Obj) 若变量 var 是 Obj 代表的类型，给出 1，表示“真”whos 给出所有 MATLAB 内存变量的属性【例 2.1-5】数据对象及其识别指令的使用。（1）cleara=1;b=2;c=3;d=4; %
12、产生 4 个数值变量Mn=a,b;c,d % 利用已赋值变量构成数值矩阵Mc=a,b;c,d % 字符串中的 a,b,c,d 与前面输入的数值变量无关Ms=sym(Mc) % Ms 是一个符号矩阵，它与前面各变量无关（2）5SizeMn=size(Mn)SizeMc=size(Mc)SizeMs=size(Ms) （3）CMn=class(Mn)CMc=class(Mc)CMs=class(Ms) （4）isa(Mn,double)isa(Mc,char)isa(Ms,sym) （5）whos Mn Mc Ms 2.2 符号数字及表达式的操作2.2.1 数值数字与符号数字之间的转换一 数值数字
13、向符号数字的转换在符号运算中，“数值类数字”会自动地转换为符号数字。亦可借助 sym 函数：sym(Num,r)或 sym(Num) 数值类数字 Num 的广义有理表达sym(Num,d) 数值类数字 Num 的 “十进制浮点”近似表达sym(Num,e) 数值类数字 Num 的带 eps 误差的理性近似表达sym(Num,f) 数值类数字 Num 的 “十六进制浮点”近似表达sym(Num)和 sym(Num)区别问题IEEE-754 标准的浮点数格式二 符号数字向双精度数字转换double(Num_sym) 把符号数字 Num_sym 转换为双精度数字2.2.2 符号数字的任意精度计算di
14、gits 显示当前环境下符号数字“十进制浮点”表示的有效数字位数digits(n) 设定符号数字“十进制浮点”表示的有效数字位数（默认 32 位）xs=vpa(x) 据表达式 x 得到 digits 指定精度下的符号数字 xsxs=vpa(x,n) 据表达式 x 得到 n 位有效数字的符号数字 xs【例 2.2-1】digits, vpa 指令的使用。digitsp0=sym(1+sqrt(5)/2)6pr=sym(1+sqrt(5)/2)%pd=sym(1+sqrt(5)/2,d)%e32r=vpa(abs(p0-pr) )e16=vpa(abs(p0-pd),16) %计算误差e32d=v
15、pa(abs(p0-pd)Digits = 32p0 =(1+sqrt(5)/2pr =7286977268806824*2(-52)pd =1.6180339887498949025257388711907e32r =.543211520368251e-16e16 =0.e32d =.543211520368251e-16 2.2.3 符号表达式的基本操作collect（合并同类项）factor（进行因式分解）numden（提取公因式）等最常用：simple(EXPR) 把 EXPR 转换成最简形式【例 2.2-2】简化 。32816xfsyms xf=(1/x3+6/x2+12/x+8)(
16、1/3);g1=simple(f)g2=simple(g1)syms xf=(1/x3+6/x2+12/x+8)(1/3);g1=simple(f)g2=simple(g1) g1 =7(2*x+1)/xg2 =2+1/x 2.2.4 表达式中的置换操作一 子表达式置换操作RS,ssub=subexpr(S,ssub) 运用符号变量 ssub 置换子表达式，并重写 S 为 RS【例 2.2-3】对符号矩阵 进行特征向量分解。dcbaclear allsyms a b c d WV,D=eig(a b;c d) % V：特征向量阵 D：特征值阵RVD,W=subexpr(V;D,W) %对矩阵元
17、素中的公共子表达式进行置换表达 二 通用置换指令RES=subs(ES,old,new) 用 new 置换 ES 中的 old 后产生 RESRES=subs(ES,new) 用 new 置换 ES 中的自由变量后产生 RES【例 2.2-4】用简单算例演示 subs 的置换规则。（1）产生符号函数syms a x;f=a*sin(x)+5 f =a*sin(x) + 5 （2）符号表达式置换f1=subs(f,sin(x),sym(y) %class(f1) （3）符号常数置换f2=subs(f,a,x,2,sym(pi/3) % a 被双精度数字置换， x 被符号数字置换class(f2)
18、 （4）双精度数值置换f3=subs(f,a,x,2,pi/3) %class(f3) （5）数值数组置换之一f4=subs(subs(f,a,2),x,0:pi/6:pi) %class(f4) （6）数值数组置换之二f5=subs(f,a,x,0:6,0:pi/6:pi) %class(f5) 82.3 符号微积分2.3.1 极限和导数的符号计算大学本科高等数学中的大多数微积分问题，都能用符号计算解决，手工笔算演绎的烦劳都可以由计算机完成。limit(f,v,a) 求极限 )(limvfa【例 2.3-1】试求 。kxx1lisyms x kLim_f=limit(1-1/x)(k*x),
19、x,inf) diff(f,v,n) 求 (n 缺省时，默认 n=1)ndvf【例 2.3-2】求 求 ， ， 。xttaflcos3df2tfdxfsyms a t x;f=a,t3;t*cos(x), log(x);df=diff(f) f 对 x 的导数dfdt2=diff(f,t,2) f 对 x 的二阶导数dfdxdt=diff(diff(f,x),t) 二阶混合导数【例 2.3-3】求 的 Jacobian（雅可比）矩阵 。)sin(co),(21212xexxf 23121xfxfsyms x1 x2;f=x1*exp(x2);x2;cos(x1)*sin(x2);v=x1 x2
20、;fjac=jacobian(f,v) 2.3.2 序列/级数的符号求和symsum(f,v,a,b) 求 f 在变量 v 取遍a, b中所有整数时的和。a,b 缺省时默认求和区间0, v-91。【例 2.3-8】求 ， 。103,ttk12)1(,(kksyms k t;f1=t k3;f2=1/(2*k-1)2,(-1)k/k;s1=simple(symsum(f1) % f1 的自变量被确认为 ts2=simple(symsum(f2,1,inf) % f1 的自变量被确认为 ksyms x y;f1=y x3;f2=1/(2*y-1)2,(-1)y/y;s1=simple(symsum
21、(f1) % f1 的自变量被确认为 xs2=simple(symsum(f2,1,inf) % f1 的自变量被确认为 y 2.3.3 符号积分int(f,v) 求 f 对变量 v 的不定积分int(f,v,a,b) 求 f 对变量 v 的定积分【例 2.3-9】求 。dx1clearsyms xf=sqrt(1+x)/x)/xs=int(f,x)s=simple(simple(s) f =(x + 1)/x)(1/2)/xs =- 2*(1/x + 1)(1/2) - 2*atan(1/x + 1)(1/2)*i)*is =- 2*(1/x + 1)(1/2) - 2*atan(1/x +
22、 1)(1/2)*i)*i 【例 2.3-10】求 。dxxbasin12syms a b x;f=a*x,b*x2;1/x,sin(x);disp(The integral of f is);pretty(int(f) 【例 2.3-11】求积分 。2 1 22)(xydzyxsyms x y zF2=int(int(int(x2+y2+z2,z,sqrt(x*y),x2*y),y,sqrt(x),x2),x,1,2)VF2=vpa(F2) %积分结果用 32 位数字表示【例 2.3-12】求阿基米德（Archimedes）螺线 在 到 间的曲线长)0(,ar度函数，并求出 时的曲线长度。2
23、,1a10（1）syms a r theta1 phi1 positivex=r*cos(theta1);x=subs(x,r,a*theta1);y=r*sin(theta1);y=subs(y,r,a*theta1);dLdth=sqrt(diff(x,theta1)2+diff(y,theta1)2);L=simple(int(dLdth,theta1,0,phi1) （2）L_2pi=subs(L,a,phi1,sym(1,2*pi)L_2pi_vpa=vpa(L_2pi) （3）L1=subs(L,a,sym(1);ezplot(L1*cos(phi1),L1*sin(phi1),0
24、,2*pi)grid onhold onx1=subs(x,a,sym(1);y1=subs(y,a,sym(1);h1=ezplot(x1,y1,0,2*pi);set(h1,Color,r,LineWidth,5)title( )legend(螺线长度 -幅角曲线, 阿基米德螺线)hold off -5 0 5 10 15 20-16-14-12-10-8-6-4-2024xy基基基基 -基基基基基基基基基基111说出以下三条指令产生的结果各属于哪种数据类型，是“双精度”对象，还是“符号”对象？3/7+0.1, sym(3/7+0.1), vpa(sym(3/7+0.1)答案c1 =0.5
25、286c2 =37/70c3 =.528571428571428571428571428571432在不加专门指定的情况下，以下符号表达式中的哪一个变量被认为是独立自由变量。sym(sin(w*t) , sym(a*exp(-X) ) , sym(z*exp(j*theta)答案ans =w ans =a ans =z 123求符号矩阵 的行列式值和逆，所得结果应采用“子表达式置换”32311aA简洁化。答案A = a11, a12, a13 a21, a22, a23 a31, a32, a33DA =a11*a22*a33-a11*a23*a32-a21*a12*a33+a21*a13*a
26、32+a31*a12*a23-a31*a13*a22IAs =(a22*a33-a23*a32)/d, -(a12*a33-a13*a32)/d, -(-a12*a23+a13*a22)/d-(a21*a33-a23*a31)/d, (a11*a33-a13*a31)/d, -(a11*a23-a13*a21)/d(-a21*a32+a22*a31)/d, -(a11*a32-a12*a31)/d, (a11*a22-a12*a21)/dd =a11*a22*a33-a11*a23*a32-a21*a12*a33+a21*a13*a32+a31*a12*a23-a31*a13*a22 4对函数
27、 ，当 为正实数时，求 。（实际上，这0)(kakfka0)(kkzf就是根据定义求 Z 变换问题。）答案Z1 =-z/(-z+a)5对于 ，求 。（提示：理论结果为 ）0x1202kkxxln答案13s_ss =log(x) 6（1）通过符号计算求 的导数 。（2）然后根据此结果，求 和ttysin)(dty 0tdy。2tdy答案d =abs(1,sin(t)*cos(t)d0_ =-1dpi_2 =0 7求出 的具有 64 位有效数字的积分值。dxexsin7.15答案matlab 2008asi =1.087849499412904913166671875948174520895458
28、535212845987519414167 8计算二重积分 。212)(xdyx答案r =1006/105 9在 区间，画出 曲线，并计算 。2,0dtxyx0sin)( )5.4(y答案y =sinint(x)y5 =1.6541404143792439835039224868515140 1 2 3 4 5 600.20.40.60.811.21.41.61.8xsinint(x)10 求 的一般积分表达式，并计算 的 32 位有效数字表xdnyn20si)( )31(y达。答案yn =1/2*pi(1/2)*gamma(1/2+1/2*n)/gamma(1+1/2*n)11 有序列 ，
29、，（在此 ， ），求这两个序列的卷积kax)(kbh)(0ba。knny0答案y1 =(-(b/a)(k+1)+1)*ak*a/(-b+a) 或y2 =(-b*bk+ak*a)/(-b+a)12 设系统的冲激响应为 ，求该系统在输入 ， 作用下的输teh3)(ttucos)(0出。（提示：运用卷积进行计算）答案hut =-3/10/exp(t)3+3/10*cos(t)+1/10*sin(t) 1513 求 的 Fourier 变换。0,)(tAetf答案F =2*A*a/(a2+w2)14 求 的 Fourier 变换，并画出 时的幅频谱。ttAtf01)( 2,A答案Fws =4*A/ta
30、o/w2*sin(1/2*tao*w)2Fw2 =4/w2*sin(w)2-6 -4 -2 0 2 4 600.511.522.533.54w4 abs(1/w2 sin(w)2)15 求 的 Laplace 反变换。463)(2ssF答案f =1/3*exp(-t)*(-2*cos(3(1/2)*t)+3(1/2)*sin(3(1/2)*t)+2) 16 利用符号运算证明 Laplace 变换的时域求导性质： 。)0()(ftLsdtf16答案Ldy =s*laplace(f(t),t,s)-f(0) 17 求 的 Z 变换表达式。Tkef )(答案F_z =z*exp(-lambda*T)
31、/(z-exp(-lambda*T)2 18 求方程 的解。2,12xyx答案x =-1/2*(1/2*5(1/2)+1/2*i*3(1/2)3+1/4*5(1/2)+1/4*i*3(1/2)-1/2*(1/2*5(1/2)-1/2*i*3(1/2)3+1/4*5(1/2)-1/4*i*3(1/2)-1/2*(-1/2*5(1/2)+1/2*i*3(1/2)3-1/4*5(1/2)+1/4*i*3(1/2)-1/2*(-1/2*5(1/2)-1/2*i*3(1/2)3-1/4*5(1/2)-1/4*i*3(1/2)y =1/2*5(1/2)+1/2*i*3(1/2)1/2*5(1/2)-1/2
32、*i*3(1/2)-1/2*5(1/2)+1/2*i*3(1/2)-1/2*5(1/2)-1/2*i*3(1/2) 19 求图 p2-1 所示信号流图的系统传递函数，并对照胡寿松主编“自动控制原理”中的例2-21 结果，进行局部性验证。图 p2-1答案传递函数 Y2U 为G1 G4 (G3 G2 + G3 G2 H5 + G3 G6 G5 + G7 G5)/(G1 G2 G4 H3 G3 H5+ G1 G4 H2 G3 H4 G5 + G1 G4 H3 G3 G6 G5 + G1 G2 G4 H3 G3+ G1 G4 H3 G7 G5 - G4 H2 H1 G7 G5 + G4 H2 G3 H
33、5 + G1 H4 G5 + G2 H1 H5+ G4 H2 G3 + H1 G6 G5 + G2 H1 + H5 + 1)局部性验证用的传递函数17a f (e + e d + c b)-e g d + g b c + e g + d + 120 采用代数状态方程法求图 p2-2 所示结构框图的传递函数 和 。UYW图 p2-2答案传递函数 Y2U 为G1 G2-G1 G2 H2 + G1 G2 H1 + 1传递函数 Y2W 为-G1 G2 H1 - G2 G3 - 1-G1 G2 H2 + G1 G2 H1 + 121 求微分方程 的通解，并绘制任意常数为 1 时解的图形。045xy答案y =1/2*(-5*x2+4*C1)(1/2)-1/2*(-5*x2+4*C1)(1/2)yy =1/2*(-5*x2+4)(1/2)-1/2*(-5*x2+4)(1/2)18-2 -1 0 1 2-2-1.5-1-0.500.511.52xy(1)y(2)22 求一阶微分方程 的解。2)(,2btax答案x =1/3*a*t3+1/2*b*t2+23 21/3 a t + 1/2 b t + 2 23 求边值问题 的解。1)0(,)(,4,3 gfgfdxgfdx答案f =exp(3*t)*sin(4*t)g =exp(3*t)*cos(4*t)
