4.6符号运算

Matlab提供了符号数学工具包（Symbolic Math Toolbox），它是操作和解决符号表达式的符号数学工具箱（函数）集合，符号表达式的运算有复合、简化、微分、积分以及求解代数方程和微分方程等。求解线性代数的问题，如求解逆、行列式的精确结果、符号矩阵的特征值。工具箱还支持可变精度运算，即支持符号计算并能以指定的精度返回结果。

在数值计算默认环境中，任何输入的数字都被映射为浮点体系中数加以保存，只能保证有限位有效数字的精度上与用户期望的精准数字相一致。Matlab符号计算环境中，保证产生的符号数字与用户期望的精准数字完全一致。符号运算的时间较长，而数值型运算速度快。

4.6.1符号运算基础

1.创建符号变量、常量和表达式

Matlab提供了sym、syms命令来定义符号变量和表达式。其常用格式：

var = sym('var')

syms var1 var2 ... varN

例如：

>> a=sym('a');

>> b=sym('beta');

>> syms a b c d x;                                          % b被重新定义

>> sym('2*sqrt(5)-1')

ans =

2*5^(1/2) - 1

>> e=[a b c; b c a; c a b]

e =

[ a, b, c]

[ b, c, a]

[ c, a, b]

>> f=a*x^3+b*x^2+c*x+d

f =

a*x^3 + b*x^2 + c*x + d

>> whos                                                         %列出当前工作空间中的变量及信息

 Name     Size           Bytes Class   Attributes

 a         1x1             112 sym               

 ans         1x1               112 sym               

 b         1x1             112 sym               

 c         1x1            112 sym               

 d         1x1            112 sym               

 e          3x3            112 sym               

 f          1x1             112 sym               

 x        1x1             112 sym     

说明：例中定义了符号变量a、b，这与执行“a='a';b='beta';”显示的值相同，但数据类型不同，前者为sym object，而后者为char array，在工作空间可以查看到数据类型不同。

2.基本运算

由于Matlab采用了重载技术，使得符号表达式的运算符和基本函数都与数值计算的几乎完全相同，使符号运算的编程很方便。

（1）符号运算中的运算符

符号基本运算符中，“+”、“-”、“*”、“\”、“/”、“^”分别实现矩阵的加、减、乘、左除、右除和求幂运算。“.*”、“.\”、“./”、“.^”分别实现符号数组（针对每个元素）的左除、右除和求幂运算。“'”、“.'”分别实现符号矩阵的共轭转置和非共轭转置。

在符号对象中，没有“大于”、“大于等于”、“小于”、“小于等于”等概念，只有“等于”、“不等于”的概念。“==”、“~=”分别用于对运算符两边的符号对象进行比较。当为“真”时，结果用1表示，当为“假”时，结果用0表示。

（2）函数

有关三角函数、双曲函数和反三角函数的符号运算与数值计算使用方法相同。指数、对数函数sqrt、exp、log2、log10等都可用于符号运算。矩阵函数diag、triu、tril、rank、eig等都可用于符号计算，但函数svd稍有区别。

3.变换函数

下面介绍几个变换函数。

sym可以将一个数值标量或矩阵转换为符号形式。其格式为：

S = sym(A,flag)

说明：sym函数中的参数flag有四种选项可以是'f'、'r'、'e'或'd'，默认值是'r'。例如：

>> a=hilb(3)

a =

   1.0000   0.5000   0.3333

   0.5000   0.3333   0.2500

   0.3333   0.2500   0.2000

>> sym(a)

ans =

[ 1, 1/2, 1/3]

[ 1/2, 1/3, 1/4]

[ 1/3, 1/4, 1/5]

>> sym(3*pi/4,'f')                                 %将该变量改为浮点数输出

ans =

2652839157010665/1125899906842624

>> sym(3*pi/4,'r')                              %将该变量改为有理数形式输出

ans =

(3*pi)/4

>> sym(3*pi/4,'e')                              %有理数形式输出，误差估计形式

ans =

(3*pi)/4 - (103*eps)/249

>> sym(3*pi/4,'d')                             %十进制小数输出，有效位数由digits定义

ans =

2.3561944901923448369984726014081

double将符号矩阵可任意精度表示的矩阵转换成双精度矩阵。例如：

>> gold=sym('(1+sqrt(5))/2')

gold =

5^(1/2)/2 + 1/2

>> double(gold)

ans =

   1.6180

eval可以计算符号表达式的值。例如：

>> syms f x; f=sym('1+x^2')

f =1+x^2

>> x=sym('2')       

>> p=eval(f)

p =

    5

sym2poly可将符号多项式变换成一般多项式，系数向量形式。poly2sym将一般多项式转换成符号表达式。例如：

>> syms f x; f=sym('2*x^2+x^3-3*x+5');

>> fp=sym2poly(f)

fp =

    1    2   -3    5

>> a=[1 2 3 4 5];

>> poly2sym(a)

ans =

x^4 + 2*x^3 + 3*x^2 + 4*x + 5

>> poly2sym(a,'s')

ans=

      s^3+2*s^2-3*s+5

4.符号计算结果的绘图

将符号计算的结果绘制成图形有两个方法：一种是根据获得的符号表达式或符号数值，转换成数值数据，再利用数值绘图函数绘制成图形。例如：

syms x y;

fx=2-3/(1+exp(x));         

xk=0:0.1:5;                

fxk=subs(fx,x,xk);                                   %获得xk函数值，结果依是符号变量

plot(xk,fxk) ;                          

另一种是利用符号绘图函数直接绘图。数值绘图函数与符号绘图函数相比，对绘图对象的操控能力更强，但符号绘图函数更简单方便。常用的绘图函数有ezplot、ezplot3、ezsurf。例如：

ezplot(fx,[0 8]);

5.其他常用符号运算函数

下面简单介绍一组关于符号运算的常用函数。见表4-4。

表4-4符号运算常用函数列表

有关表4-4中符号运算函数的使用，举例如下：

>> syms x y ;

>> f1=collect(x^3-2*x^2+1+3*x^3+3*x^2)

f1 =

4*x^3 + x^2 + 1

>> f2=factor(sym('x^4-1'))                              %

f2 =

[ x - 1, x + 1, x^2 + 1]

>> f3=simplify(exp(log(x+y)))                         %化简

f3 =

x + y

>> [n,d]=numden(x/y+y/x)                           %通分，结果为

n =

x^2 + y^2

d =

x*y

4.6.2线性代数中的符号运算

1．基本代数运算

基本代数运算举例如下：

>> syms t ;a=sym([cos(2*t),sin(3*t);-sin(t),cos(4*t)])

a =

[ cos(2*t), sin(3*t)]

[ -sin(t), cos(4*t)]

>> syms a11 a12 a21 a22

>> b1=sym([1 2;3 4]);b2=[a11,a12;a21,a22];

>> b1*b2

ans =

[  a11 + 2*a21,  a12 + 2*a22]

[ 3*a11 + 4*a21, 3*a12 + 4*a22]

>> c=sym(magic(3))

c =

[ 8, 1, 6]

[ 3, 5, 7]

[ 4, 9, 2]

>> rank(c)

ans =

3

>> trace(c)

ans =

15

2.特征值和特征向量的计算

求取矩阵的特征多项式使用charpoly函数，求特征值、特征向量用eig函数。例如：

>> h=sym( [3 -1 1;2 0 1 ;1 -1 2])

h =

[ 3, -1, 1]

[ 2,      0, 1]

[ 1, -1, 2]

>> charpoly(h)                                   %返回特征多项式系数

ans =

[ 1, -5, 8, -4]

>> [v,d]=eig(h)                                       %求特征值与特征向量

v=

[ 0, 1]

[ 1, 1]

[ 1, 0]

d =

[ 1, 0, 0]

[ 0, 2, 0]

[ 0, 0, 2]

4.6.3符号微分与积分运算

1.符号变量的确定

当符号表达式中含有多于一个的变量时，通常解题时，将围绕自由符号变量进行。如果不指定，Matlab将按默认规则次序选择，大写字母比小写字母靠后，Matlab将按照与小写字母x的ASCII码距离自动识别自由符号变量，小写字母i、j不能成为自由符号变量。可利用函数findsym查询Matlab在符号表达式中哪一个变量被认为是自由符号变量。其格式是：

findsym(S,N)

说明：S输入符号表达式，N用以设定查找与之相关的N个变量。使用方法举例如下：

>> syms a b t u v w x y z;

>> s=3*a*y+b*z;

>> findsym(s)                                        %查找函数f中全部的符号变量      

ans =

a,b,y,z

>> findsym(s,1)                                    %查找函数f中自由符号变量

ans =

y

>> f=u+v+w+x+y+z;

>> findsym(f,5)                                      %查找函数f中排在前五的符号变量

ans =

x,y,w,z,v

2．微分运算

diff函数可用来求导符号表达式的微分。使用时可以指定符号表示式、符号变量及求导的阶数。diff命令也可以使用符号矩阵作为参数，微分对矩阵的每个元素进行，求导的符号变量应相同。例如：

>> p=sym('a*x^2+b*x+c')

p =

a*x^2 + b*x + c

>> diff(p)                                             % p对x求导

ans =

b + 2*a*x

>> diff(p,'a')                                          % p对a求导

ans =

x^2

>> diff(p,2)                                           % p对x求二阶导数

ans =

2*a

>> k=sym([3/2,sin(2*x+1);4/x^2,exp(3*x)])

>> diff(k)

ans =

[     0,     2*cos(2*x + 1)]

[ -8/x^3,    3*exp(3*x)   ]

3．极限运算

当某一极限存在时，可利用函数limit求取极限。其格式：

limit (F,x,a)

limit (F,x,a,'right')

limit (F,x,a, 'left')

说明：对x求趋近a时的极限，省略a，表示求趋近0的极限，'right'表示从右趋近于a，'left'表示从左趋近于a。

[例4-30]求当n→inf,求limit(1+1/n)^1/2的值。

>> syms n ;

>> limit((1+1/n)^(1/2),n,inf)

ans =

1

再举例说明使用方法。

>> syms n x;

>> limit((1+x/n)^n,n,inf)

ans =

exp(x)

>> limit(sin(x)/x,x,0)

ans =

1

>> limit(1/x,x,0)

ans =

NaN

>> limit(1/x,x,0,'left'), limit(1/x,x,0,'right')

ans =

-Inf

ans =

Inf

4．积分运算

函数的积分运算包括求函数的不定积分和定积分，积分比微分复杂得多。积分或逆求导不一定是以封闭形式存在，或许存在但软件也许找不到，或者软件可求解，但超过内存或时间限制，当Matlab不能找到逆导数时，它将返回未经计算的命令。其基本格式：

int(f,x)

int(f,x,a,b)

说明：表示对符号变量x求从a到b的定积分。缺省x，表示对自由符号变量求不定积分，a，b是数值也可是是符号变量。如果f是矩阵，如同函数diff一样，积分函数int对符号数组的每一个元素进行运算。例如：

>> syms x s m n; f=sin(s+2*x);      

>> int(f)                                              %以x为积分变量进行积分

ans=

      -cos(s + 2*x)/2

>> int(f,s)                                           %以s为积分变量进行积分

ans =

      -cos(s+2*x)

>> int(f,pi/2,pi)                              %以x为积分变量从/2到进行积分

ans=

      -cos(s)

>> int(f,'m','n')                                   %求x从m到n的积分

ans =

      cos(2*m + s)/2 - cos(2*n + s)/2

>> syms x; int(exp(-1*x^2),-inf,inf)

ans =

pi^(1/2)

5．求积运算和taylor级数展开

（1）求积运算

当符号变量的级数和存在时，使用symsum函数来求和。symsum即可以对常数级数求和，也可以对含变量的级数求和。函数求表达式的级数和较通用格式如下：

symsum(f,x,a,b)

说明：计算表达式f的级数和，x为自变量，x省略表示对默认自由符号变量求和，a、b为x的取值范围。

解：由题意知求解。

>>syms n; symsum( 1/(2*n-1)^2,1,inf)

ans =

pi^2/8

若求此级数前五项的和，即求

>> syms n;symsum( 1/(2*n-1)^2,1,5)

ans =

117469/99225

[例4-32]求含变量的级数和。

>> syms k x;symsum( x^k,'k',0,inf)

ans =

piecewise([1 <= x, Inf], [abs(x) < 1, -1/(x - 1)])

（2）taylor级数展开

Matlab中提供了将函数展开为幂级数的函数taylor，其较通用格式如下：

taylor(f,x,a)

说明：f为符号表达式，x为自变量，对f进行taylor级数展开，参数a指定函数f在自变量x=a处展开，a的默认值为0。

>> syms x;f=taylor((1-3*x+x^2)^(1/3))

f =

- (31*x^5)/9 - (16*x^4)/9 - x^3 - (2*x^2)/3 - x + 1

>> syms x; f=log(x+1);g=cos(x);

>> taylor(f,x,1)

ans =

x/2 + log(2) - (x - 1)^2/8 + (x - 1)^3/24 - (x - 1)^4/64 + (x - 1)^5/160 - 1/2

>> taylor(g,x,'order',10)

ans =

x^8/40320 - x^6/720 + x^4/24 - x^2/2 + 1

4.6.4符号方程和微分方程求解

1.简单代数方程

代数方程是指未涉及微积分运算的方程，在Matlab中，求解用符号表达式表示的代数方程可由函数solve实现。如果表达式不是一个方程式(不含等号)，则在求解之前函数solve将表达式置成等于0。函数格式如下：

y = solve(eqn,var)

说明：求解符号表达式表示的代数方程eqn，求解变量为var。求解代数方程组格式如下：

[y1,...,yN] = solve(eqns,vars)

说明：表示求解符号表达式组成的方程组eqns，求解多个变量vars，若不指定求解变量，则由默认规则确定。举例如下：

>>syms a b c x; solve( ' a*x^2+b*x+c ' )          

ans =

-(b + (b^2 - 4*a*c)^(1/2))/(2*a)

-(b - (b^2 - 4*a*c)^(1/2))/(2*a)

结果是符号向量，其元素是方程的2个解。如果想对非缺省x变量求解，solve必须指定变量。

>> f=solve('a*x^2+b*x+c=0','a')

f =

-(c + b*x)/x^2

>> solve( ' cos(x)=sin(x) ' )                  

ans=

      1/4*pi

求解方程组。例如：

>> [x,y]=solve( 'x+3*y=7','x^2+y=3' )

x =

   1

-2/3

y =

   2

23/9

>> a=solve( 'x+3*y=7','x^2+y=3' )             %a是结构体变量

a =

   x: [2x1 sym]

   y: [2x1 sym]

>> a.x

ans =

   1

-2/3

2.微分方程和微分方程组

常微分方程有时很难求解，Matlab提供了功能强大的工具，函数dsolve计算常微分方程的符号解。因为我们要求解微分方程，就需要用一种方法将微分包含在表达式中。所以，dsolve语法与大多数其它函数有一些不同，用大写字母D来表示求微分，D2、D3等表示重复求微分，并以此来设定方程。任何D后所跟的字母为因变量。方程用符号表达式D2y=0来表示。独立变量可以指定或由symvar规则选定。函数dsolve格式如下：

S = dsolve(eqn,cond)

[y1,...,yN] = dsolve(eqns,conds)

Y = dsolve(eqns,conds)

>> dsolve( ' Dy=1+y^2 ' )

ans =

tan(C6 + t)

       i

         -i

其中，C6是任意常数。若不给出初值条件，则求方程的通解，若加入初值y(0)=1条件，则可求得特解。

>> dsolve(' Dy=1+y^2 ',' y(0)=1 ')

ans =

tan(pi/4 + t)

>> y=dsolve(' D2y=cos(2*x)-y ',' Dy(0)=0 ',' y(0)=1 ')

y=

cos(2*x) - cos(t)*(cos(2*x) - 1)

[例4-36]应用dsolve求解例4-27微分方程，并验证dsolve得到的解析结果与数值求值的ode45解算结果是一致的。

>> dsolve('D2y+y=1-1/(2*pi)*t^2','y(0)=0,Dy(0)=0')

ans =

1/pi - cos(t)*(1/pi + 1) - t^2/(2*pi) + 1

编写脚本文件：

t0=0;tf=3*pi;

[t,y]=ode45('odefun1',[t0,tf], [0;0]);                 %函数odefun1在例4-27中建立完成

ys=1/pi - cos(t)*(1/pi + 1) - t.^2/(2*pi) + 1;      %注意与dsolve给出解析表达式书写的不同

plot(t,ys,t,y(:,1),'o');

legend('解析解','数值积分解');

                           图4-12 数值积分解与解析解

结果如图4-12所示，可以看到dsolve得到的解析结果与的ode45解算结果是一致的。