ch05符号计算.doc

第五章 符号计算

符号计算的特点：一，运算以推理解析的方式进行，因此不受计算误差积累问题困扰；二，符号计算，或给出完全正确的封闭解，或给出任意精度的数值解（当封闭解不存在时）；三，符号计算指令的调用比较简单，经典教科书公式相近；四，计算所需时间较长，有时难以忍受。
在MATLAB中，符号计算虽以数值计算的补充身份出现，但涉及符号计算的指令使用、运算符操作、计算结果可视化、程序编制以及在线帮助系统都是十分完整、便捷的。
MATLAB的升级和符号计算内核Maple的升级，决定着符号计算工具包的升级。但从用户使用角度看，这些升级所引起的变化相当细微。即使这样，本章还是及时作了相应的更新和说明。如MATLAB 6.5+ 版开始启用Maple VIII的计算引擎，从而克服了Maple V计算“广义Fourier变换”时的错误（详见第5.4.1节）。
5.1 符号对象和符号表达式
5.1.1 符号对象的生成和使用
【例5.1.1-1】符号常数形成中的差异
a1=[1/3,pi/7,sqrt(5),pi+sqrt(5)] % <1>
a2=sym([1/3,pi/7,sqrt(5),pi+sqrt(5)]) % <2>
a3=sym([1/3,pi/7,sqrt(5),pi+sqrt(5)],'e') % <3>
a4=sym('[1/3,pi/7,sqrt(5),pi+sqrt(5)]') % <4>
a24=a2-a4
a1 =
0.3333 0.4488 2.2361 5.3777
a2 =
[ 1/3, pi/7, sqrt(5), 6054707603575008*2^(-50)]
a3 =
[ 1/3-eps/12, pi/7-13*eps/165, sqrt(5)+137*eps/280, 6054707603575008*2^(-50)]
a4 =
[ 1/3, pi/7, sqrt(5), pi+sqrt(5)]
a24 =
[ 0, 0,
0, 189209612611719/35184372088832-pi-5^(1/2)]

【例5.1.1-2】演示：几种输入下产生矩阵的异同。
a1=sym([1/3,0.2+sqrt(2),pi]) % <1>
a2=sym('[1/3,0.2+sqrt(2),pi]') % <2>
a3=sym('[1/3 0.2+sqrt(2) pi]') % <3>
a1_a2=a1-a2 %
a1 =
[ 1/3, 7269771597999872*2^(-52), pi]
a2 =
[ 1/3, 0.2+sqrt(2), pi]
a3 =
[ 1/3, 0.2+sqrt(2), pi]
a1_a2 =
[ 0, 1.4142135623730951010657008737326-2^(1/2), 0]

【例5.1.1-3】把字符表达式转换为符号变量
y=sym('2*sin(x)*cos(x)')
y=simple(y)
y =
2*sin(x)*cos(x)
y =
sin(2*x)

【例5.1.1-4】用符号计算验证三角等式。
syms fai1 fai2;y=simple(sin(fai1)*cos(fai2)-cos(fai1)*sin(fai2))
y =
sin(fai1-fai2)

【例5.1.1-5】求矩阵的行列式值、逆和特征根
syms a11 a12 a21 a22;A=[a11,a12;a21,a22]
DA=det(A),IA=inv(A),EA=eig(A)
A =
[ a11, a12]
[ a21, a22]
DA =
a11*a22-a12*a21
IA =
[ a22/(a11*a22-a12*a21), -a12/(a11*a22-a12*a21)]
[ -a21/(a11*a22-a12*a21), a11/(a11*a22-a12*a21)]
EA =
[ 1/2*a11+1/2a22+1/2*(a11^2-2*a11*a22+a22^2+4*a12*a21)^(1/2)]
[ 1/2*a11+1/2*a22-1/2*(a11^2-2*a11*a22+a22^2+4*a12*a21)^(1/2)]

【例5.1.1-6】验证积分。
syms A t tao w;yf=int(A*exp(-i*w*t),t,-tao/2,tao/2);Yf=simple(yf)
Yf =
2*A*sin(1/2*tao*w)/w

5.1.2 符号计算中的算符和基本函数
5.1.3 识别对象类别的指令
【例5.1.3-1】数据对象及其识别指令的使用。
（1）
clear,a=1;b=2;c=3;d=4;
Mn=[a,b;c,d]
Mc='[a,b;c,d]'
Ms=sym(Mc)
Mn =
1 2
3 4
Mc =
[a,b;c,d]
Ms =
[ a, b]
[ c, d]

（2）
SizeMn=size(Mn),SizeMc=size(Mc),SizeMs=size(Ms)
SizeMn =
2 2
SizeMc =
1 9
SizeMs =
2 2

（3）
CMn=class(Mn),CMc=class(Mc),CMs=class(Ms)
CMn =
double
CMc =
char
CMs =
sym

（4）
isa(Mn,'double'),isa(Mc,'char'),isa(Ms,'sym')
ans =
1
ans =
1
ans =
1

（5）
whos Mn Mc Ms
Name Size Bytes Class

Mc 1x9 18 char array
Mn 2x2 32 double array
Ms 2x2 408 sym object

Grand total is 21 elements using 458 bytes
5.1.4 符号表达式中自由变量的确定
【例5.1.4-1】对独立自由符号变量的自动辨认。
（1）
syms a b x X Y;k=sym('3');z=sym('c*sqrt(delta)+y*sin(theta)');
EXPR=a*z*X+(b*x^2+k)*Y;

（2）
findsym(EXPR)
ans =
X, Y, a, b, c, delta, theta, x, y

（3）
findsym(EXPR,1)
ans =
x

（4）
findsym(EXPR,2),findsym(EXPR,3)
ans =
x,y
ans =
x,y,theta

【例5.1.4-2】findsym确定自由变量是对整个矩阵进行的。
syms a b t u v x y;A=[a+b*x,sin(t)+u;x*exp(-t),log(y)+v]
findsym(A,1)
A =
[ a+b*x, sin(t)+u]
[ x*exp(-t), log(y)+v]
ans =
x

5.2 符号表达式和符号函数的操作
5.2.1 符号表达式的操作
【例5.2.1-1】按不同的方式合并同幂项。
EXPR=sym('(x^2+x*exp(-t)+1)*(x+exp(-t))');
expr1=collect(EXPR)
expr2=collect(EXPR,'exp(-t)')
expr1 =
x^3+2*exp(-t)*x^2+(1+exp(-t)^2)*x+exp(-t)
expr2 =
x*exp(-t)^2+(2*x^2+1)*exp(-t)+x^2+1)*x

【例5.2.1-2】factor指令的使用
（1）
syms a x;f1=x^4-5*x^3+5*x^2+5*x-6;factor(f1)
ans =
(x-1)*(x-2)*(x-3)*(x+1)

（2）
f2=x^2-a^2;factor(f2)
ans =
(x-a)*(x+a)

（3）
factor(1025)
ans =
5 5 41

【例5.2.1-3】对多项式进行嵌套型分解
clear;syms a x;f1=x^4-5*x^3+5*x^2+5*x-6;horner(f1)
ans =
-6+(5+(5+(-5+x)*x)*x)*x

【例5.2.1-4】写出矩阵各元素的分子、分母多项式
（1）
syms x;A=[3/2,(x^2+3)/(2*x-1)+3*x/(x-1);4/x^2,3*x+4];
[n,d]=numden(A)
pretty(simplify(A)) %<3>
n =
[ 3, x^3+5*x^2-3]
[ 4, 3*x+4]
d =
[ 2, (2*x-1)*(x-1)]
[ x^2, 1]

[ 3 2 ]
[ x + 5 x - 3 ]
[3/2 -----------------]
[ (2 x - 1) (x - 1)]
[ ]
[ 4 ]
[---- 3 x + 4 ]
[ 2 ]
[ x ]

（2）
pretty(simplify(n./d))

【例5.2.1-5】简化
（1）
syms x;f=(1/x^3+6/x^2+12/x+8)^(1/3);
sfy1=simplify(f),sfy2=simplify(sfy1)
sfy1 =
((2*x+1)^3/x^3)^(1/3)
sfy2 =
((2*x+1)^3/x^3)^(1/3)

（2）
g1=simple(f),g2=simple(g1)
g1 =
(2*x+1)/x
g2 =
2+1/x

【例5.2.1-6】简化
syms x;ff=cos(x)+sqrt(-sin(x)^2);
ssfy1=simplify(ff),ssfy2=simplify(ssfy1)
ssfy1 =
cos(x)+(sin(x)^2)^(1/2)
ssfy2 =
cos(x)+(-sin(x)^2)^(1/2)

gg1=simple(ff),gg2=simple(gg1)
gg1 =
cos(x)+i*sin(x)
gg2 =
exp(i*x)

5.2.2 符号函数的求反和复合
【例5.2.2-1】求的反函数
syms x;f=x^2;g=finverse(f)
Warning: finverse(x^2) is not unique.
> In D:\MATLAB6P5\toolbox\symbolic\@sym\finverse.m at line 43
g =
x^(1/2)
fg=simple(compose(g,f)) %验算g(f(x))是否等于x
fg =
x

【例5.2.2-2】求的复合函数
（1）
syms x y u fai t;f=x/(1+u^2);g=cos(y+fai);fg1=compose(f,g)
fg1 =
cos(y+fai)/(1+u^2)

（2）
fg2=compose(f,g,u,fai,t)
fg2 =
x/(cos(y+t)^2+1)

5.2.3 置换及其应用
5.2.3.1 自动执行的子表达式置换指令
【例5.2.3.1-1】演示子表达式的置换表示。
clear all,syms a b c d W;[V,D]=eig([a b;c d]);
[RVD,W]=subexpr([V;D],W) %<2>
RVD =
[ -(1/2*d-1/2*a-1/2*W)/c, -(1/2*d-1/2*a+1/2*W)/c]
[ 1, 1]
[ 1/2*d+1/2*a+1/2*W, 0]
[ 0, 1/2*d+1/2*a-1/2*W]
W =
(d^2-2*a*d+a^2+4*b*c)^(1/2)

5.2.3.2 通用置换指令
【例5.2.3.2-1】用简单算例演示subs的置换规则。
（1）
syms a x;f=a*sin(x)+5;

（2）
f1=subs(f,'sin(x)',sym('y')) %<2>
f1 =
a*y+5

（3）
f2=subs(f,{a,x},{2,sym(pi/3)}) %<3>
f2 =
3^(1/2)+5

（4）
f3=subs(f,{a,x},{2,pi/3}) %<4>
f3 =
6.7321

（5）
f4=subs(subs(f,a,2),x,0:pi/6:pi) %<5>
f4 =
5.0000 6.0000 6.7321 7.0000 6.7321 6.0000 5.0000

（6）
f5=subs(f,{a,x},{0:6,0:pi/6:pi}) %<6>
f5 =
5.0000 5.5000 6.7321 8.0000 8.4641 7.5000 5.0000

5.2.4 符号数值精度控制和任意精度计算
5.2.4.1 向双精度数值转换的doblue指令
5.2.4.2 任意精度的符号数值
【例5.2.4.2-1】指令使用演示。
digits
Digits = 32

p0=sym('(1+sqrt(5))/2');
p1=sym((1+sqrt(5))/2)
e01=vpa(abs(p0-p1))
p1 =
7286977268806824*2^(-52)
e01 =
.543211520368250e-16

p2=vpa(p0)
e02=vpa(abs(p0-p2),64)
p2 =
1.6180339887498948482045868343657
e02 =
.61882279690820194237137864551377e-31

digits
Digits = 32

5.2.5 符号对象与其它数据对象间的转换
【例5.2.5-1】符号、数值间的转换。
phi=sym((1+sqrt(5))/2)
double(phi)
phi =
7286977268806824*2^(-52)
ans =
1.6180

【例5.2.5-2】各种多项式表示形式之间的转换
syms x;f=x^3+2*x^2-3*x+5;
sy2p=sym2poly(f)
p2st=poly2str(sy2p,'x')
p2sy=poly2sym(sy2p)
pretty(f,'x')
sy2p =
1 2 -3 5
p2st =
x^3 + 2 x^2 - 3 x + 5
p2sy =
x^3+2*x^2-3*x+5

3 2
x + 2 x - 3 x + 5

5.3 符号微积分
5.3.1 符号序列的求和
【例5.3.1-1】求，
syms k t;f1=[t k^3];f2=[1/(2*k-1)^2,(-1)^k/k];
s1=simple(symsum(f1))
s2=simple(symsum(f2,1,inf))
s1 =
[ 1/2*t*(t-1), k^3*t]
s2 =
[ 1/8*pi^2, -log(2)]

5.3.2 符号微分和矩阵
【例5.3.2-1】求、和
syms a t x;f=[a,t^3;t*cos(x), log(x)];
df=diff(f)
dfdt2=diff(f,t,2)
dfdxdt=diff(diff(f,x),t)
df =
[ 0, 0]
[ -t*sin(x), 1/x]
dfdt2 =
[ 0, 6*t]
[ 0, 0]
dfdxdt =
[ 0, 0]
[ -sin(x), 0]

【例5.3.2-2】求的矩阵。
syms x1 x2 x3;f=[x1*exp(x2);x2;cos(x1)*sin(x2)];
v=[x1 x2];fjac=jacobian(f,v)
fjac =
[ exp(x2), x1*exp(x2)]
[ 0, 1]
[ -sin(x1)*sin(x2), cos(x1)*cos(x2)]

5.3.3 符号积分
5.3.3.1 通用积分指令
5.3.3.2 交互式近似积分指令
5.3.3.3 符号积分示例
【例5.3.3.3-1】求。演示：积分指令对符号函数矩阵的作用。
syms a b x;f=[a*x,b*x^2;1/x,sin(x)]
disp('The integral of f is');pretty(int(f))
The integral of f is

[ 2 3]
[1/2 a x 1/3 b x ]
[ ]
[ log(x) -cos(x) ]

【例5.3.3.3-2求。演示如何使用mfun指令获取一组积分值。
（1）
F1=int('1/log(t)','t',0,'x')
F1 =
-Ei(1,-log(x))

（2）
x=0.5:0.1:0.9
F115=-mfun('Ei',1,-log(x))
x =
0.5000 0.6000 0.7000 0.8000 0.9000
F115 =
-0.3787 -0.5469 -0.7809 -1.1340 -1.7758

【例5.3.3.3-3】求积分。注意：内积分上下限都是函数。
syms x y z
F2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2)
VF2=vpa(F2)
F2 =
1610027357/6563700-6072064/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2^(3/4)
VF2 =
224.92153573331143159790710032805

【例5.3.3.3-4】利用rsums求积分。（与例5.3.3.3-2结果比较）
syms x positive;px=0.5/log(0.5*x);rsums(px)

图5.3-1

5.3.4 符号卷积
【例5.3.4-1】本例演示卷积的时域积分法：已知系统冲激响应，求输入下的输出响应。
syms T t tao;ut=exp(-t);
ht=exp(-t/T)/T;
uh_tao=subs(ut,t,tao)*subs(ht,t,t-tao);
yt=int(uh_tao,tao,0,t);
yt=simple(yt)
yt =
-(exp(-t)-exp(-t/T))/(T-1)

【例5.3.4-2】本例演示通过变换和反变换求取卷积。系统冲激响应、输入同上例，求输出。
对式(53.4-1)两边进行Laplace变换得，因此有
syms s;yt=ilaplace(laplace(ut,t,s)*laplace(ht,t,s),s,t);yt=simple(yt)
yt =
(-exp(-t)+exp(-t/T))/(T-1)

【例5.3.4-3】求函数和的卷积。
syms tao;t=sym('t','positive');
ut=sym('Heaviside(t)-Heaviside(t-1)');ht=t*exp(-t);
yt=int(subs(ut,t,tao)*subs(ht,t,t-tao),tao,0,t);
yt=collect(yt,'Heaviside(t-1)')
yt =
(-1+exp(1-t)*t)*Heaviside(t-1)+1+(-t-1)*exp(-t)

5.4 符号积分变换
5.4.1 Fourier变换及其反变换
【例 5.4.1-1】求的Fourier变换。本例演示三个重要内容：单位阶跃函数和单位脉冲函数的符号表示；fourier指令的使用；simple指令的表现。
（1）求Fourier变换
syms t w;ut=sym('Heaviside(t)'); % <1>
UT=fourier(ut)
UTC=maple('convert',UT,'piecewise','w') % <3>
UTS=simple(UT)
UT =
pi*Dirac(w)-i/w
UTC =
PIECEWISE([undefined, w = 0],[0, otherwise])
UTS =
-i/w

（2）求Fourier反变换进行验算
Ut=ifourier(UT,w,t)
Uts=ifourier(UTS,w,t)
Ut =
1/2+1/2*Heaviside(t)-1/2*Heaviside(-t)
Uts =
1/2*Heaviside(t)-1/2*Heaviside(-t)

【例5.41-2】用fourier指令求例5.1.1-6中方波脉冲的Fourier变换。本例演示：fourier, simple 指令的配合使用。
（1）
syms A t w
syms tao positive % <2>
yt=sym('Heaviside(t+tao/2)-Heaviside(t-tao/2)');
Yw=fourier(A*yt,t,w)
Ywc=maple('convert',Yw,'piecewise','w') %计算结果起指示作用 <5>
Yws=simple(Yw)
Yw =
A*(exp(1/2*i*tao*w)*(pi*Dirac(w)-i/w)-exp(-1/2*i*tao*w)*(pi*Dirac(w)-i/w))
Ywc =
-i*A*exp(1/2*i*tao*w)/w+i*A*exp(-1/2*i*tao*w)/w
Yws =
2*A*sin(1/2*tao*w)/w

（2）
Yt=ifourier(Yw,w,t)
Yst=ifourier(Yws,w,t)
Yt =
1/2*A*(Heaviside(t+1/2*tao)-Heaviside(-t-1/2*tao)-Heaviside(t-1/2*tao)+Heaviside(-t+1/2*tao))
Yst
1/2*A*(Heaviside(t+1/2*tao)-Heaviside(-t-1/2*tao)-Heaviside(t-1/2*tao)+Heaviside(-t+1/2*tao))

【例5.4.1-3】求的Fourier变换，在此是参数，是时间变量。本例演示：fourier的缺省调用格式的使用要十分谨慎；在被变换函数中包含多个符号变量的情况下，对被变换的自变量给予指明，可保证计算结果的正确。
syms t x w;ft=exp(-(t-x))*sym('Heaviside(t-x)');
F1=simple(fourier(ft,t,w))
F2=simple(fourier(ft))
F3=simple(fourier(ft,t))
F1 =
1/exp(i*x*w)/(1+i*w)
F2 =
i*exp(-i*t*w)/(i+w)
F3 =
i*exp(-t*(2+i*t))/(i+t)

5.4.2 Laplace变换及其反变换
【例5.4.2-1】求的Laplace变换。
syms t s;syms a b positive %<1>
Dt=sym('Dirac(t-a)'); %<2>
Ut=sym('Heaviside(t-b)'); %<3>
Mt=[Dt,Ut;exp(-a*t)*sin(b*t),t^2*exp(-t)];MS=laplace(Mt,t,s)
MS =
[ exp(-a*s), exp(-b*s)/s]
[ b/((s+a)^2+b^2), 2/(1+s)^3]

【例5.4.2-2】验证Laplace时移性质： 。
syms t s;t0=sym('t0','positive'); %<1>
ft=sym('f(t-t0)')*sym('Heaviside(t-t0)')
FS=laplace(ft,t,s),FS_t=ilaplace(FS,s,t)
ft =
f(t-t0)*Heaviside(t-t0)
FS =
exp(-s*t0)*laplace(f(t),t,s)
FS_t =
f(t-t0)*Heaviside(t-t0)

5.4.3 Z变换及其反变换
【例5.4.3-1】求序列 的Z变换，并用反变换验算。
（1）
syms n
Delta=sym('charfcn[0](n)'); % <2>
D0=subs(Delta,n,0); %
D15=subs(Delta,n,15); %
disp('[D0,D15]');disp([D0,D15])
[D0,D15]
[ 1, 0]

（2）求序列的Z变换
syms z;fn=2*Delta+6*(1-(1/2)^n);FZ=simple(ztrans(fn,n,z));
disp('FZ = ');pretty(FZ),FZ_n=iztrans(FZ,z,n)
FZ =
2
4 z + 2
--------------
2
2 z - 3 z + 1
FZ_n =
2*charfcn[0](n)+6-6*(1/2)^n

5.5 符号代数方程的求解
5.5.1 线性方程组的符号解
【例 5.5.1-1】求线性方程组的解。本例演示，符号线性方程组的基本解法。
该方程组的矩阵形式是 。该式简记为。求符号解的指令如下
A=sym([1 1/2 1/2 -1;1 1 -1 1;1 -1/4 -1 1;-8 -1 1 1]);
b=sym([0;10;0;1]);X1=A\b
X1 =
[ 1]
[ 8]
[ 8]
[ 9]

【例 5.5.1-2】求解上例前3 个方程所构成的“欠定”方程组，并解释解的含义。
A2=A(1:3,:);X2=A2\b(1:3,1)
syms k;
XX2=X2+k*null(A2)
A2*XX2
X2 =
[ 0]
[ 8]
[ 4]
[ 6]
XX2 =
[ 1/3*k]
[ 8]
[ 4+4/3*k]
[ 6+k]
ans =
[ 0]
[ 10]
[ 0]

5.5.2 一般代数方程组的解
【例5.5.2-1】求方程组，关于的解。
S=solve('u*y^2+v*z+w=0','y+z+w=0','y','z')
disp('S.y'),disp(S.y),disp('S.z'),disp(S.z)
S =
y: [2x1 sym]
z: [2x1 sym]
S.y
[ -1/2/u*(-2*u*w-v+(4*u*w*v+v^2-4*u*w)^(1/2))-w]
[ -1/2/u*(-2*u*w-v-(4*u*w*v+v^2-4*u*w)^(1/2))-w]
S.z
[ 1/2/u*(-2*u*w-v+(4*u*w*v+v^2-4*u*w)^(1/2))]
[ 1/2/u*(-2*u*wv-(4*u*w*v+v^2-4*u*w)^(1/2))]

【例5.5.2-2】用solve指令重做例 5.5.1-2。即求，，构成的“欠定”方程组解。
syms d n p q;eq1=d+n/2+p/2-q;eq2=n+d+q-p-10;eq3=q+d-n/4-p;
S=solve(eq1,eq2,eq3,d,n,p,q);S.d,S.n,S.p,S.q
Warning: 3 equations in 4 variables.
> In D:\MATLAB6P5\toolbox\symbolic\solve.m at line 110
In D:\MATLAB6P5\toolbox\symbolic\@sym\solve.m at line 49
ans =
d
ans =
8
ans =
4*d+4
ans =
3*d+6

【例5.5.2-3】求的解。
clear all,syms x;s=solve('(x+2)^x=2','x')
s =
.69829942170241042826920133106081

5.6 符号微分方程的求解
5.6.1 符号解法和数值解法的互补作用
5.6.2 求微分方程符号解的一般指令
5.6.3 微分方程符号解示例
【例5.6.3-1】求的解。
S=dsolve('Dx=y,Dy=-x');
disp([blanks(12),'x',blanks(21),'y']),disp([S.x,S.y])
x y
[ cos(t)*C1+sin(t)*C2, -sin(t)*C1+cos(t)*C2]

【例5.6.3-2】图示微分方程的通解和奇解的关系。
y=dsolve('y=x*Dy-(Dy)^2','x')
clf,hold on,ezplot(y(2),[-6,6,-4,8],1)
cc=get(gca,'Children'); % <3>
set(cc,'Color','r','LineWidth',5) % <4>
for k=-2:0.5:2;ezplot(subs(y(1),'C1',k),[-6,6,-4,8],1);end
hold off,title('\fontname{隶书}\fontsize{16}通解和奇解')
y =
[ x*C1-C1^2]
[ 1/4*x^2]

图 5.6-1

【例5.6.3-3】求解两点边值问题： 。（注意：相应的数值解法比较复杂）。
y=dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=