实验MATLAB符号计算.docx
《实验MATLAB符号计算.docx》由会员分享，可在线阅读，更多相关《实验MATLAB符号计算.docx（33页珍藏版）》请在冰点文库上搜索。
实验MATLAB符号计算
实验四符号计算
符号计算的特点：
一，运算以推理解析的方式进行，因此不受计算误差积累问题困扰；二，符号计算，或给出完全正确的封闭解，或给出任意精度的数值解（当封闭解不存在时）；三，符号计算指令的调用比较简单，经典教科书公式相近；四，计算所需时间较长，有时难以忍受。
在MATLAB中，符号计算虽以数值计算的补充身份出现，但涉及符号计算的指令使用、运算符操作、计算结果可视化、程序编制以及在线帮助系统都是十分完整、便捷的。
MATLAB的升级和符号计算内核Maple的升级，决定着符号计算工具包的升级。
但从用户使用角度看，这些升级所引起的变化相当细微。
即使这样，本章还是及时作了相应的更新和说明。
如MATLAB6.5+版开始启用MapleVIII的计算引擎，从而克服了MapleV计算“广义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.33330.44882.23615.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/30.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】用符号计算验证三角等式
。
symsfai1fai2;y=simple（sin（fai1）*cos（fai2）-cos（fai1）*sin（fai2））
y=
sin（fai1-fai2）
【例5.1.1-5】求矩阵
的行列式值、逆和特征根
symsa11a12a21a22;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/2*a22+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】验证积分
。
symsAttaow;yf=int（A*exp（-i*w*t）,t,-tao/2,tao/2）;Yf=simple（yf）
Yf=
2*A*sin（1/2*w*tao）/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=
12
34
Mc=
[a,b;c,d]
Ms=
[a,b]
[c,d]
（2）
SizeMn=size（Mn）,SizeMc=size（Mc）,SizeMs=size（Ms）
SizeMn=
22
SizeMc=
19
SizeMs=
22
（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）
whosMnMcMs
NameSizeBytesClass
Mc1x918chararray
Mn2x232doublearray
Ms2x2312symobject
Grandtotalis21elementsusing362bytes
5.1.4符号表达式中自由变量的确定
【例5.1.4-1】对独立自由符号变量的自动辨认。
（1）
symsabxXY;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确定自由变量是对整个矩阵进行的。
symsabtuvxy;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）
symsax;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=
-（a-x）*（a+x）
（3）
factor（1025）
ans=
5541
【例5.2.1-3】对多项式进行嵌套型分解
clear;symsax;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）
symsx;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]
[32]
[x+5x-3]
[3/2-----------------]
[（2x-1）（x-1）]
[]
[4]
[----3x+4]
[2]
[x]
（2）
pretty（simplify（n./d））
[32]
[x+5x-3]
[3/2-----------------]
[（2x-1）（x-1）]
[]
[4]
[----3x+4]
[2]
[x]
【例5.2.1-5】简化
（1）
symsx;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】简化
symsx;ff=cos（x）+sqrt（-sin（x）^2）;
ssfy1=simplify（ff）,ssfy2=simplify（ssfy1）
ssfy1=
cos（x）+（-1+cos（x）^2）^（1/2）
ssfy2=
cos（x）+（-1+cos（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】求
的反函数
symsx;f=x^2;g=finverse（f）
g=
x^（1/2）
fg=simple（compose（g,f））%验算g（f（x））是否等于x
fg=
x
【例5.2.2-2】求
的复合函数
（1）
symsxyufait;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/（1+cos（y+t）^2）
5.2.3置换及其应用
5.2.3.1自动执行的子表达式置换指令
【例5.2.3.1-1】演示子表达式的置换表示。
clearall,symsabcdW;[V,D]=eig（[ab;cd]）;
[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）
symsax;f=a*sin（x）+5;
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.00006.00006.73217.00006.73216.00005.0000
（6）
f5=subs（f,{a,x},{0:
6,0:
pi/6:
pi}）%<6>
f5=
5.00005.50006.73218.00008.46417.50005.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=
.543211520368251e-16
p2=vpa（p0）
e02=vpa（abs（p0-p2）,64）
p2=
1.6180339887498948482045868343656
e02=
.38117720309179805762862135448622e-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】各种多项式表示形式之间的转换
symsx;f=x^3+2*x^2-3*x+5;
sy2p=sym2poly（f）
p2st=poly2str（sy2p,'x'）
p2sy=poly2sym（sy2p）
pretty（f,'x'）
sy2p=
12-35
p2st=
x^3+2x^2-3x+5
p2sy=
x^3+2*x^2-3*x+5
5.3符号微积分
5.3.1符号序列的求和
【例5.3.1-1】求
，
symskt;f1=[tk^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】求
、
和
symsatx;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】求
的
矩阵。
symsx1x2x3;f=[x1*exp（x2）;x2;cos（x1）*sin（x2）];
v=[x1x2];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】求
。
演示：
积分指令对符号函数矩阵的作用。
symsabx;f=[a*x,b*x^2;1/x,sin（x）];
disp（'Theintegraloffis'）;pretty（int（f））
Theintegraloffis
[23]
[1/2ax1/3bx]
[]
[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.50000.60000.70000.80000.9000
F115=
-0.3787-0.5469-0.7809-1.1340-1.7758
【例5.3.3.3-3】求积分
。
注意：
内积分上下限都是函数。
symsxyz
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=
64/225*2^（3/4）-6072064/348075*2^（1/2）+14912/4641*2^（1/4）+1610027357/6563700
VF2=
224.92153573331143159790710032805
【例5.3.3.3-4】利用rsums求
积分。
（与例5.3.3.3-2结果比较）
symsxpositive;px=0.5/log（0.5*x）;rsums（px）
图5.3-1
5.3.4符号卷积
【例5.3.4-1】本例演示卷积的时域积分法：
已知系统冲激响应
，求
输入下的输出响应。
symsTttao;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=
-1/（T-1）/exp（t）+1/（T-1）/exp（t/T）
【例5.3.4-2】本例演示通过变换和反变换求取卷积。
系统冲激响应、输入同上例，求输出。
对式（5.3.4-1）两边进行Laplace变换得
，因此有
symss;yt=ilaplace（laplace（ut,t,s）*laplace（ht,t,s）,s,t）;yt=simple（yt）
yt=
（exp（-t/T）-exp（-t））/（T-1）
【例5.3.4-3】求函数
和
的卷积。
symstao;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=
（exp（1-t）*t-1）*heaviside（t-1）+1+（-t-1）/exp（t）
5.4符号积分变换
5.4.1Fourier变换及其反变换
【例5.4.1-1】求
的Fourier变换。
本例演示三个重要内容：
单位阶跃函数和单位脉冲函数的符号表示；fourier指令的使用；simple指令的表现。
（1）求Fourier变换
symstw;ut=sym（'Heaviside（t）'）;%<1>
UT=fourier（ut）
UTC=maple（'convert',UT,'piecewise','w'）%<3>
UTS=simple（UT）
UT=
2*pi*dirac（w）
UTC=
PIECEWISE（[pi*NaN,w=0],[0,otherwise]）
UTS=
2*pi*dirac（w）
（2）求Fourier反变换进行验算
Ut=ifourier（UT,w,t）
Uts=ifourier（UTS,w,t）
Ut=
1
Uts=
1
【例5.4.1-2】用fourier指令求例5.1.1-6中方波脉冲的Fourier变换。
本例演示：
fourier,simple指令的配合使用。
（1）
symsAtw
symstaopositive%<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*（i*exp（-1/2*i*tao*w）/w+pi*dirac（w））
Ywc=
PIECEWISE（[A*（i*exp（-1/2*i*tao*w）/w+pi*NaN）,w=0],[i*A*exp（-1/2*i*tao*w）/w,otherwise]）
Yws=
A*（i*exp（-1/2*i*tao*w）/w+pi*dirac（w））
（2）
Yt=ifourier（Yw,w,t）
Yst=ifourier（Yws,w,t）
Yt=
A*heaviside（-t+1/2*tao）
Yst=
A*heaviside（-t+1/2*tao）
【例5.4.1-3】求
的Fourier变换，在此
是参数，
是时间变量。
本例演示：
fourier的缺省调用格式的使用要十分谨慎；在被变换函数中包含多个符号变量的情况下，对被变换的自变量给予指明，可保证计算结果的正确。
symstxw;ft=exp（-（t-x））*sym（'Heaviside（t-x）'）;
F1=simple（fouri
升级会员 