第22章　符号数学工具 MATLAB所具有的符号数学工具箱与其它所有工具不同，它适用于广泛的用途，而不是针对一些特殊专业或专业分支。另外，MATLAB符号数学工具箱与其它的工具箱区别还因为它使用字符串来进行符号分析，而不是基于数组的数值分析。为此，本章包含了该工具箱的教学辅导材料。 22.1 引言 符号数学工具箱是操作和解决符号表达式的符号数学工具箱(函数)集合，有复合、简化、微分、积分以及求解代数方程和微分方程的工具。另外还有一些用于线性代数的工具，求解逆、行列式、正则型式的精确结果，找出符号矩阵的特征值而无由数值计算引入的误差。工具箱还支持可变精度运算，即支持符号计算并能以指定的精度返回结果。 符号数学工具箱中的工具是建立在功能强大的称作Maple软件的基础上。它最初是由加拿大的滑铁卢（Waterloo）大学开发的。当要求MATLAB进行符号运算时，它就请求Maple去计算并将结果返回到MATLAB命令窗口。因此，在MATLAB中的符号运算是MATLAB处理数字的自然扩展。 22.2 符号表达式 符号表达式是代表数字、函数、算子和变量的MATLAB字符串，或字符串数组。不要求变量有预先确定的值，符号方程式是含有等号的符号表达式。符号算术是使用已知的规则和给定符号恒等式求解这些符号方程的实践，它与代数和微积分所学到的求解方法完全一样。符号矩阵是数组，其元素是符号表达式。 MATLAB在内部把符号表达式表示成字符串，以与数字变量或运算相区别；否则，这些符号表达式几乎完全象基本的MATLAB命令。表22.1列有几则符号表达式例子以及MATLAB等效表达式。 表22.1 符号表达式 MATLAB表达式    ' 1/(2*x^n) '   y= ' 1/sqrt(2*x) '    ' cos(x^2)-sin(2*x) '   M=sym( ' [a，b；c，d] ' )   f=int( ' x^3/sqrt(1-x) ' ， ' a ' ， ' b ' )   MATLAB符号函数可让用户用多种方法来操作这些表达式，比如， >> diff( ' cos(x) ' ) % differentiate cos(x) with respect to x ans= -sin(x) >> M=sym( ' [a，b；c，d] ' ) % create a symbolic matrix M M= [a，b] [c，d] >> determ(M) % find the determinant of the symbolic matrix M ans= a*d-b*c 请注意，上面的第一个例子的符号表达式是用单引号以隐含方式定义的。它告诉MATLAB ' cos(x) ' 是一个字符串并说明diff( ' cosx ' )是一个符号表达式而不是数字表达式；然而在第二个例子中，用函数sym显式地告诉MATLAB M=sym( ' [a，b；c，d] ' )是一符号表达式。在MATLAB可以自己确定变量类型的场合下，通常不要求显式函数sym。 正如在第八章所阐述，MATLAB中函数function argument形式是与function( ' argument ' )等价的。其中，function是一个函数，argument是一字符串。例如，MATLAB可以构造diff cos(x)和diff( ' cos(x) ' )两者都意味diff (sym ' cos(x) ' )。但第一种形式显然更便于输入。然而，很多时候sym是必要的。在上述的第二个例子中， >> M=[a，b；c，d] % M is a numeric matrix using value of a through d ???Uundefine function or variable a. >> M= ' [a，b；c，d] ' % M is a character string， but not a symbolic matrix M= [a，b；c，d] >> M=sym( ' [a，b；c，d] ' ) % M is a symbolic matrix M= [a，b] [c，d] M以三种方式定义： 数字型(如果a、b、c、d已预先确定)、字符串型或符号矩阵型。许多符号函数非常巧妙能够自动将字符转变为符号表达式。但在某些情况下，尤其是建立符号数组时，必须用函数sym，特别地将字符串变为符号表达式。隐含形式，例如diff cos(x)，对于那些不需要参考先前结果的简单任务，最有用。但是最简单形式(无引号)要求一个参量，它是一个单字符的字符串、不包含插入的空格。 >> diff x^2+3*x+5 % the argument is equivalent to ' x^2+3*x+5 ' ans= 2*x+3 >> diff x^2 + 3*x + 5 % spaces break the argument into separate strings ???Error using==>diff Too manyinput arguments 无变量的符号表达式称作符号常量。符号常量常常与整数很难区别，例如 >> f=symop( ' (3*4-2)/5+1 ' ) % reduce a symbolic constant to its simplest form f = 3 >> isstr(f) % is f a string? (1=yes， 0=no) ans= 1 在这个例子中，f代表符号常数 ' 3 ' ；而不是数字3。正如第六章所阐述的，MATLAB是以字符ASCII码形式来存储字符串的。所以，如果对字符串进行数字运算，则在运算中，采用各字符串的ASCII码值。因为数字51是字符 ' 3 ' 的ASCII表示，所以f加1在数值上不能得到期望的结果 >> f+1 ans= 52 符号变量 当字符表达式中含有多于一个的变量时，只有一个变量是独立变量。如果不告诉MATLAB哪一个变量是独立变量，MATLAB将基于以下规则选择一个： 在符号表达式中缺省的独立变量是唯一的，除去i和j的小写字母，不是单词的一部分。如果没有这种字母，就选择x作为独立变量。如字符不是唯一的，就选择在字母顺序中最接近x的字母。如果有相连的字母，就选择在字母表中较后的那一个。 缺省的独立变量，有时称作自由变量，在表达式 ' 1/(5+cos(x)) ' 中是 ' x ' ；在 ' 3*y+z ' 中是 ' y ' ；在 ' a+sin(t) ' 是 ' t ' 。在表式 ' sin(pi/4)-cos(3/5) ' 中自由符号变量是 ' x ' ，因为此式是一个符号常数无符号变量。可利用函数symvar询问MATLAB在符号表达式中哪一个变量它认为是独立变量。 >> symvar( ' a*x+y*) % find the default symbolic variable ans= x >> symvar( ' a*t+s/(u+3) ' ) % u is the closest to ' x ' ans= u >> symvar( ' sin(omega) ' ) % ' omega ' is not a singlee character。 ans= x >> symvar( ' 3*i+4*j ' ) % i and j are equel to sqrt(-1) ans= x >> symvar( ' y+3*s ' ， ' t ' ) % find the variable closest to t rather than x ans= s 如果利用规则symvar不能找到一个缺省独立变量，它便假定无独立变量并返回x。这一结论对含有由多个字母组成的变量，如：alpha或s2的表达式，或不含变量的符号常数均成立。如果需要，绝大多数命令都使用用户选项以指定独立变量。 >> diff( ' x^n ' ) % differentiate with respect to the default variable ' x ' ans= x^n*n/x >> diff( ' x^n ' ， ' n ' ) % differentiate x^n with respect to ' n ' ans= x^n*log(x) >> diff( ' sin(omega) ' ) % differentiate using the default variables (x) ans= 0 >> diff( ' sin(omega) ' ， ' omega ' ) % specify the independent variable ans= cos(omega) 22.3 符号表达式运算 一旦创建了一个符号表达式，或许想以某些方式改变它；也许希望提取表达式的一部分，合并两个表达式或求得表达的数值。有许多符号工具可以帮助完成这些任务。 所有符号函数(很少特殊例外的情况，讨论于后)作用到符号表达式和符号数组，并返回符号表达式或数组。其结果有时可能看起来象一个数字，但事实上它是一个内部用字符串表示的一个符号表达式。正如我们前面所讨论的，可以运用MATLAB函数isstr来找出像似数字的表达式是否真是一个整数或是一个字符串。 提取分子和分母 如果表达式是一个有理分式(两个多项式之比)，或是可以展开为有理分式(包括哪些分母为1的分式)，可利用numden来提取分子或分母。例如，给定如下的表达式：  在必要时，numden将表达式合并、有理化并返回所得的分子和分母。进行这项运算的MATLAB语句是： >> m= ' x^2 ' % create a simple expression m= x^2 >> [n，d]=numden(m) % extract the numerator and denominator n= x^2 d= 1 >> f= ' a*x^2/(b-x) ' % create a rational expression f= a*x^2/(b-x) >> [n，d]=numden(f) % extract the numerator and denominator n= a*x^2 d= b-x 前二个表达式得到期望结果。 >> g= ' 3/2*x^2+2/3*x-3/5 ' % rationalize and extract the parts g= 3/2*x^2+2/3*x-3/5 >> [n，d]=numden(g) n= 45*x^2+20*x-18 d= 30 >> h= ' (x^2+3)/(2*x-1)+3*x/(x-1) ' % the sum of rational polynomials h= (x^2+3)/(2*x-1)+3*x/(x-1) >> [n，d]=numden(h) % rationalize and extract n= x^3+5*x^2-3 d= (2*x-1)*(x-1) 在提取各部分之前，这二个表达式g和h被有理化，并变换成具有分子和分母的一个简单表达式。 >> k=sym( ' [3/2，(2*x+1)/3；4/x^2，3*x+4] ' ) % try a symbolic array k= [ 3/2，(2*x+1)/3] [4/x^2， 3*x+4] >> [n，d]=numden(k) n= [3， 2*x+1] [4， 3*x+4] d= [ 2，3] [x^2，1] 这个表达式k是符号数组，numden返回两个新数组n和d，其中n是分子数组，d是分母数组。如果采用s=numden（f)形式，numden仅把分子返回到变量s中。 标准代数运算 很多标准的代数运算可以在符号表达式上执行，函数symadd、symsub、symlnul和symdiv为加、减、乘、除两个表达式，sympow将一个表达式上升为另一个表达式的幂次。例如： 给定两个函数  >> f= ' 2*x^2+3*x-5 ' % define the symbolic expression f= 2*x^2+3*x-5 >> g= ' x^2-x+7 ' g= x^2-x+7 >> symadd(f，g) % find an expression for f+g ans= 3*x^2+2*x+2 >> symsub(f，g) % find an expression for f-g ans= x^2+4*x-12 >> symmul(f，g) % find an expression for f*g ans= (2*x^2+3*x-5)*(x^2-x+7) >> symdiv(f，g) % find an expression for f/g ans= (2*x^2+3*x-5)/(x^2-x+7) >> sympow(f， ' 3*x ' ) % find an expression for  ans= (2*x^2+3*x-5)^3** 另一个通用函数可让用户用其它的符号变量、表达式和算子创建新的表达式。symop取由逗号隔开的、多至16个参量。各个参量可为符号表达式、数值或算子(' + '、' - '、'*'、' / '、' ^ '、' ( '或' ) ')，然后symop可将参量联接起来，返回最后所得的表达式. >> f= ' cos(x) ' % create an expression f= cos(x) >> g= ' sin(2*x) ' % create another expression g= sin(2*x) >> symop(f，'/ ＇，g，'+'，3) % combine them ans= cos(x)/sin(2*x)+3 所有这些运算也同样用数组参量进行。 高级运算 MATLAB具有对符号表达式执行更高级运算的功能。函数compose把f(x)和g(x)复合成f(g(x))。函数finverse求表达式的函数逆，而函数symsum求表达式的符号和。 给定表达式  >> f= ' 1/(1+x^2) ' ； % create the four expression >> g= ' sin(x) ' ； >> h= ' 1/(1+u^2) ' ； >> k=' sin(v) ' ； >> compose(f，g) % find an expression for f(g(x)) ans= 1/(1+sin(x)^2) >> compose(g，f) % find an expression for g(f(x)) ans= sin(1/(1+x^2)) compose也可用于含有不同独立变量的函数表达式。 >> compose(h，k，'u'，'v') % given h(u)，k(v)，find(k(v)) ans= 1/(1+sin(v)^2) 表达式譬如f(x)的函数逆g(x)，满足g(f(x))=x。例如，的函数逆是ln(x)，因为ln()=x。sin(x)的函数逆是arcsin(x)，函数的函数逆是arcsin。函数fincerse返回表达式的函数逆。如果解不是唯一就给出警告。 >> finverse( ' 1/x) % the inverse of 1/x is 1/x since ' 1/（1/x)=x ' ans= 1/x >> finverse( ' x^2 ' ) % g(x^2)=x has more than one solution Warning： finverse(x^2) is not unique ans= x^(1/2) >> finverse( ' a*x+b ' ) % find the solution to ' g(f(x))=x ' ans= -(b-x)/a >> finverse( ' a*b+c*d-a*z ' )， ' a ' ) % find the solution to ' g(f(a))=a ' ans= -(c*d-a)/(b-z) symsun函数求表达式的符号和有四种形式：symsun(f)返回；symsum(f， ' s ' )返回，symsun(f，a，b)返回；最普通的形式symsun(f， ' s ' ，a，b)返回。 让我们试一试，它应返回： 。 >> symsum('x^2') ans= 1/3*x^3-1/2* x^2+1/6*x 又怎么样呢?它应返回。 >> sym('(2*n-1)^2'，1，'n') ans= 11/3*n+8/3-4*(n+1)^2+4/3*(n+1)^3 >> factor(ans) % change the form ( we will revisit 'factor' later on) ans= 1/3*n*(2*n-1)*(2*n+1) 最后让我们试一试，其返回应是。 >> symsum( ' 1/(2*n-1)^2 ' ，1，inf) ans= 1/8*pi^2 变换函数 本节提出许多工具，将符号表达式变换成数值或反之。有极少数的符号函数可返回数值。然而请注意，某些符号函数能自动地将一个数字变换成它的符号表达式，如果该数字是函数许多参量中的一个。 函数sym可获取一个数字参量并将其转换为符号表达式。函数numneric的功能正好相反，它把一个符号常数（无变量符号表达式）变换为一个数值。 >> phi=' (1+sqrt(5))/2 ' % the ' golden ' ratio phi= (1+sqrt(5))/2 % convert to a numeric value >> numeric(phi) ans= 1.6180 正如第六章所介绍，函数eval将字符串传给MATLAB以便计算。所以eval是另一个可用于把符号常数变换为数字或计算表达式的函数。 >> eval(phi) % execute the string ' (1+sqrt(5))/2 ans= 1.6180 正如所期望那样，numeric和eval返回相同数值。 符号函数sym2poly将符号多项式变换成它的MATLAB等价系数向量。函数poly2syrn功能正好相反，并让用户指定用于所得结果表达式中的变量。 >> f=' 2*x^2+x^3-3*x+5 ' % f is the symbolic polynomials f= 2*x^2+x^3-3*x+5 >> n=sym2poly(f) % extract eht numeric coefficient vector n= 1 2 -3 5 >> poly2sym(n) % recreate the polynomials in x (the default) ans= 2*x^2+x^3-3*x+5 >> poly2sym(n，' s ') % recreate the polynomials in s ans= s^3+2*s^2-3*s+5 变量替换 假设有一个以x为变量的符号表达式，并希望将变量转换为y。MATLAB提供一个工具称作subs，以便在符号表达式中进行变量替换。其格式为subs（f，new，old)，其中f是符号表达式，new和old是字符、字符串或其它符号表达式。‘新’字符串将代替表达式f中各个‘旧’字符串。以下有几个例子： >> f= ' a*x^2+b*x+c ' % create a function f(x) f= a*x^2+b*x+c >> subs(f，' s '，' x ') % substitute ' s ' for ' x ' in the expression f ans= a*s^2+b*s+c >> subs(f，' alpha '，' a ') % substitute ' alpha ' for ' a ' in f ans= alpha*x^2+b*x+c >> g=' 3*x^2+5*x-4 ' % create another function g= 3*x^2+5*x-4 >> h=subs(g，' 2 '，' x ') % substitute ' 2 ' for ' x ' in g h= 18 >> isstr(h) % show that the result is a symbolic expression ans= 1 最后一个例子表明subs如何进行替换，并力图简化表达式。因为替换结果是一个符号常数，MATLLAB可以将其简化为一个符号值。注意，因为subs是一个符号函数，所以它返回一个符号表达式。尽管看似数字，实质上是一个符号常数。为了得到数字，我们需要使用函数numeric或eval来转换字符串。 >> numeric(h) % convert a symbolic expression to a number ans= 18 >> isstr(ans) % show that the result is a numeric value ans= 0 22.4 微分和积分 微分和积分是微积分学研究和应用的核心，并广泛地用在许多工程学科。MATLAB符号工具能帮助解决许多这类问题。 微分 符号表达式的微分以四种形式利用函数diff： >> f= ' a*x^3+x^2-b*x-c ' % define a symbolic expression f= a*x^3+x^2-b*x-c >> diff(f) % differentiate with respect to the default variable x ans= 3*a*x^2+2*x-b >> diff(f，'a ') % differentiate with respect to a ans= x^3 >> diff(f，2) % differentiate twice with respect to x ans= 6*a*x+2 >> diff(f，' a '，2) % differentiate twice with respect to a ans= 0 函数diff也可对数组进行运算。如果F是符号向量或数组，diff(F)对数组内的各个元素进行微分。 >> F=sym(' [a*x， b*x^2； c*x^3， d*s] ') % create a symbolic array F= [ a*x， b*x^2] [c*x^3， d*s] >> diff(F) % differentiate the element with respect to x ans= [ a，2*b*x] [3*c*x^2， 0] 注意函数diff也用在MATLAB，计算数值向量或矩阵的数值差分。对于一个数值向量或矩阵M，diff(M)计算M(2： m，： )-M(1： m-1，： )的数值差分，如下所示： >> m=[(1： 8).^2)] % create a vector M= 1 4 9 16 25 36 49 64 >> diff(M) % find the differences between elements ans= 3 5 7 9 11 13 15 如果diff的表达式或可变参量是数值，MATLAB就非常巧妙地计算其数值差分；如果参量是符号字符串或变量，MATLAB就对其表达式进行微分。 积分 积分函数int(f)，其中f是一符号表达式，它力图求出另一符号表达式F使diff(F)=f。正如从研究微分学所了解的，积分比微分复杂得多。积分或逆求导不一定是以封闭形式存在，或许存在但软件也许找不到，或者软件可明显地求解，但超过内存或时间限制。当MATLAB不能找到逆导数时，它将返回未经计算的命令。 >> int( ' log(x)/exp(x^2) ' ) % attempt to integrate ans= log(x)/exp(x^2) 同微分一样，积分函数有多种形式。形式int(f)相对于缺省的独立变量求逆导数；形式(f，' s ')相对于符号变量s积分；形式int(f，a，b)和int(f，' s '，a，b)，a，b是数值，求解符号表达式从a到b的定积分；形式int(f，' m ' ，' n ')和形式int(f，' s '，' m '，' n ')，其中m，n是符号变量，求解符号表达式从m到n的定积分。 >> f=' sin(s+2*x) ') % crate a symbolic function f= sin(s+2*x) >> int(f) % integrate with respect to x ans= -1/2*cos(s+2*x) >> int(f，' s ') % integrate with respect to s ans= -cos(s+2*x) >> int(f，pi/2，pi) % integrate with respect to x from /2 to ans= -cos(x) >> int(f，' s '，pi/2，pi) % integrate with respect to s from /2 to  ans= cos(2*x)-sin(2*x) >> int(f，' m '，' n ') % integrate with respect to x from m to n ans= -1/2*cos(s+2*n)+1/2*cos(s+2*m) 正如函数diff一样，积分函数int对符号数组的每一个元素进行运算。 >> F=sym( ' [a*x，b*x^2；c*x^3，d*s] ' ) % create a symbolic array F= [ a*x，b*x^2] [c*x^3， d*s] >> diff(F) % ubtegrate the array elements with respect to x ans= [1/2*a*x^2，1/3*b*x^3] [1/4*c*x^4， d*s*x] 22.5 符号表达式画图 在许多的场合，将表达式可视化是有利的。MATLAB提供了函数ezplot来完成该任务。 >> y=' 16*x^2+64*x+96 ' % expression to plot y= 16*x^2+64*x+96 >> ezplot(y)  图22.1 符号函数16*x^2+64*x+96 (-2≤x≤2)  图22.2 符号函数16*x^2+64*x+96 0≤x≤6 正如图22.1所示，ezplot绘制了定义域为-2≤x≤2的给定符号函数，并相应地调整了y轴比例，还加了网格栅和标志。在这个例子中，我们感兴趣的时间是从0到6。让我们再试一下，并指定时间范围，见图22.2 >> ezplot(y，[0 6]) % plot y for 0≤x≤6 现在，在所感兴趣的范围内显示得略好些。一旦该图处在图形窗口内，它可以象其它图象一样作修改。 22.6 符号表达式简化和格式化 有时MATLAB返回的符号表达式难以理解，有许多工具可以使表达式变得更易读懂。第一个就是函数pretty，该命令以类似于数学课本上的形式来显示符号表达式。 >> f=taylor( ' log(x+1)/(x-5) ' ) % 6 terms is the default f = -1/5*x+3/50*x^2-41/750*x^3+293/7500*x^4-1207/37500*x^5+O(x^6) >> pretty(f) 2 41 3 293 4 1207 5 6 - 1/5 x + 3/50 x - -----x + ------x - -------x + O(x ) 750 7500 37500 符号表达式可用许多等价形式来提供。在不同的场合，某种形式可能胜于另一种。MATLAB用许多命令来简化或改变符号表达式。 >> f=sym( ' (x^2-1)*(x-2)*(x-3) ' ) % create a function f= (x^2-1)*(x-2)*(x-3) >> collect(f) % collect all like terms ans= x^4-5*x^3+5*x^2+5*x-6 >> hornor(ans) % change to Horner or nested representation ans= -6+(5+(5+(-5+x)*x)*x)*x >> factor(ans) % express as a product of polynomials ans= (x-1)*(x-2)*(x-3)*(x+1) >> expand(f) % distribute products over sums ans= x^4-5*x^3+5*x^2+5*x-6 Simplify是功能强大、通用的工具。它利用各种类型代数恒等式，包括求和、积分和分数幂、三角、指数和log函数、Bessel函数、超几何函数和函数，来简化表达式。 下面例子说明函数的乘幂。 >> simplify( ' log(2*x/y) ' ) ans= log(2)+log(x)-log(y) >> simplify( ' sin(x)^2+3*x+cos(x)^2-5 ' ) ans= 3*x-4 >> simplify( ' (-a^2+1)/(1-a) ' ) ans= a+1 其中要讨论的最后一个函数是最有用的，但也是最不正统的。函数simple试用了几种不同的简化工具，然后选择在结果表达式中含有最少字符的那种形式。让我们看一下立方根：  >> f=' (1/x^3+6/x^2+12/x+8)^(1/3) ' % create the expression f= (1/x^3+6/x^2+12/x+8)^(1/3) >> simple(f) % simplify it simplify ： (2*x+1)/x ans= (2*x+1)/x >> simple(ans) % try it once again - another method may help combine(trig) 2+1/x ans= 2+1/x 正如所见，simple试用了几种可简化表达式的简化方式，并让看到每一个尝试的结果。有时，它帮助多次使用函数simple并对第一次的结果作不同的简化操作，如上所作。simple对于含有三角函数的表达式尤为有用。让我们试一下 >> simple( ' cos(x)+sqrt(-sin(x)^2) ' ) % simplify a trig expression simplify： cos(x)+(cos(x)^2-1)^(1/2) radsimp： cos(x)+i*sin(x) combine(trig)： cos(x)+(-1/2+1/2*cos(2*x))^(1/2) factor： cos(x)+(-sin(x)^2)^(1/2) expand： cos(x)+(-sin(x)^2)^(1/2) convert(exp)： 1/2*exp(i*x)+1/2/exp(i*x)+1/4*4^(1/2)*(exp(i*x)-1/exp(i*x)) convert(tan)： (1-tan(1/2*x)^2)/(1+tan(1/2*x)^2)+(-4*tan(1/2*x)^2/(1+tan(1/2*)^2)^2)^(1/2) ans = cos(x)+i*sin(x)