MATLAB符号计算.docx
《MATLAB符号计算.docx》由会员分享，可在线阅读，更多相关《MATLAB符号计算.docx（31页珍藏版）》请在冰豆网上搜索。
MATLAB符号计算
第2章符号计算
所谓符号计算是指：
解算数学表达式、方程不是在离散化的数值点上进行，而是凭借一系列恒等式，数学定理，通过推理和演绎，力求获得解析结果。
这种计算建立在数值完全准确表达和推演严格解析的基础之上，因此所得结果是完全准确的。
.1符号对象和符号表达式
.1.1基本符号对象和运算算符
101生成符号对象的基本规则
生成符号对象的MATLAB规则：
（1）必须定义
（2）包含符号对象的也一定是符号对象
102精准符号数字和符号常数
sc=sym（Num）
【例2.1-1】演示：
精准符号数字或数字表达式的创建。
1）
clearall
R1=sin（sym（0.3））%
R2=sin（sym（3e-1））%
R3=sin（sym（3/10））%
R4=sin（sym（'3/10'））%
disp（['R1属于什么类别？
答：
',class（R1）]）
disp（['R1与R4是否相等？
（是为1，否为0）答：
',int2str（logical（R1==R4））]）
R1=
sin（3/10）
R2=
sin（3/10）
R3=
sin（3/10）
R4=
sin（3/10）
R1属于什么类别？
答：
sym
R1与R4是否相等？
（是为1，否为0）答：
1
103基本符号变量
sym:
创建符合数字、符号常数 、符号表达式
syms:
创建多个符号变量、抽象符号函数以及变量
.1.2符号计算中的函数指令
符号与数值计算中有些函数同名，使用时注意数据类型！
表2.1-2MATLAB中可调用的符号计算函数指令
类别
情况描述
符号工具包函数
三角函数、双曲函数及反函数；除atan2外
指数、对数函数（如exp,expm）
复数函数（如abs,angle）
矩阵分解函数（如eig,svd）
方程求解函数solve
微积分函数（如diff,int）
积分变换和反变换函数（如laplace,ilaplace）
绘图函数（如ezplot,ezsurf）
特殊函数
单位脉冲和阶跃函数（如dirac,heaviside）
函数（如beta,gamma）
椭圆积分（如ellipke）
贝塞尔函数（如besseli,besselj）
MuPAD库函数
借助evalin和feval指令可调用比mfunlist所列范围更广泛的MuPAD库函数；需要具备MuPAD语言知识。
〖说明〗
.1.3符号表达式和符号函数
101符号表达式和符号函数
102自由符号变量
【例2.1-2】
1）
clear
symsabcxyuv%
symsF（X,Y,Z）%
k=sym（3）%
G=sym（'p*sqrt（q）+r*sin（t）'）%
EXPR=a*G*u+（b*x^2+k）*v%
f（x,y）=a*x^2+b*y^2-c%
disp（F）%
k=
3
G=
p*q^（1/2）+r*sin（t）
EXPR=
v*（b*x^2+3）+a*u*（p*q^（1/2）+r*sin（t））
f（x,y）=
a*x^2+b*y^2-c
F（X,Y,Z）
symbolicfunctioninputs:
X,Y,Z
2）
symvar（EXPR）%
%
ans=
[a,b,p,q,r,t,u,v,x]
symvar（EXPR,1）%
ans=
x
3）
disp（symvar（f））%
[a,b,c,x,y]
disp（symvar（f,2））%
[x,y]
【例2.1-3】用符号法求方程
的解。
1）
clear
symsuvwz%
Eq=u*w^2+z*w-v%
Eq=
u*w^2+z*w-v
2）
symvar（Eq,1）%<5>
ans=
w
3）
R1=solve（Eq）%<7>
R1=
-（z+（z^2+4*u*v）^（1/2））/（2*u）
-（z-（z^2+4*u*v）^（1/2））/（2*u）
4）
S1=solve（Eq,z）%<9>
S1=
（-u*w^2+v）/w
5）
disp（simplify（subs（Eq,z,S1）））%
0
〖说明〗
【例2.1-4】
symsabtuvxy
A=[a+b*x,sin（t）+u;x*exp（-t）,log（y）+v]%
symvar（A,1）%
A=
[a+b*x,u+sin（t）]
[x*exp（-t）,v+log（y）]
ans=
x
.1.4符号对象的识别
用于识别数据的指令：
class（var）给出变量var的数据类别（如double,sym等）
isa（var,'Obj'）若变量var是Obj代表的类别，给出1，表示“真”
whos给出所有MATLAB内存变量的属性
【例2.1-5】数据对象及其识别指令的使用。
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%
NameSizeBytesClassAttributes
Mc1x918char
Mn2x232double
Ms2x260sym
.1.5符号运算机理和变量假设
101符号运算的工作机理
102对符号变量的限定性假设
assume（expr,set）%set为‘integer’,‘rational’,‘real’表示整数集，有理数集，实数集
a=sym（‘a’,res）%res为‘positive’,‘real’表示正数，实数
103清除变量和撤销假设
clearx与symsxclear的区别
【例2.1-6】符号变量的默认数域是复数域。
1）
clearall%
symsx%
f=x^3+475*x/100+5/2;%
r=solve（f,x）%
r=
-1/2
（79^（1/2）*i）/4+1/4
1/4-（79^（1/2）*i）/4
assumptions（x）%
ans=
[emptysym]
2）
assume（x,'real'）%
r21=solve（f,x）%
r21=
-1/2
symsxclear%
assume（imag（x）==0）%
r22=solve（f,x）
r22=
-1/2
disp（assumptions（x））%
imag（x）==0
3）
sym（x,'clear'）%
assume（real（x）>0）%
r3=solve（f,x）
ans=
x
r3=
（79^（1/2）*i）/4+1/4
1/4-（79^（1/2）*i）/4
disp（assumptions（x））%
04）
assumeAlso（imag（x）>0）%
r4=solve（f,x）
r4=
（79^（1/2）*i）/4+1/4
disp（assumptions（x））
[0.2符号数字及表达式的操作
.2.1符号数字转换成双精度数字
Newdata=double（num_sym）
.2.2符号数字的任意精度表达形式
digits（n）设定十进制符号数字有效位数
vpa（x,n）设定精度，只在运行当时起作用
【例2.2-1】
1）
reset（symengine）%<1>
sa=sym（'1/3+sqrt
（2）'）%
sa=
2^（1/2）+1/3
2）
digits%<3>
Digits=32
formatlong
a=1/3+sqrt
（2）%
sa_Plus_a=vpa（sa+a,20）%<6>
sa_Minus_a=vpa（sa-a,20）%<7>
%
a=
1.747546895706428
sa_Plus_a=
3.4950937914128567869
sa_Minus_a=
-0.000000000000000022658064826339973669
3）
sa32=vpa（sa）%<8>
digits（48）%<9>
sa5=vpa（sa,5）%<10>
sa48=vpa（sa）%<11>
sa32=
1.747546895706428382135022057543
sa5=
1.7475
sa48=
1.74754689570642838213502205754303141190300520871
〖说明〗
.2.3符号表达式的基本操作
collect,expand,factor,horner,numden,simplify,pretty,coeffs
【例2.2-2】简化
。
symsx
f=（1/x^3+6/x^2+12/x+8）^（1/3）
g1=simplify（f）%
f=
（12/x+6/x^2+1/x^3+8）^（1/3）
g1=
（（2*x+1）^3/x^3）^（1/3）
.2.4表达式中的置换操作
101公因子法简化表达
subexpr（S,‘w’）从S中提取公因子w
102通用置换指令
subs（ES,old,new）用new置换ES中的old
【例2.2-4】
1）
clear
symsabx;
f1=a*sin（x）+b
f1=
b+a*sin（x）
2）
f2=subs（f1,sin（x）,'log（y）'）%<4>
f2=
b+a*log（y）
.3符号微积分
可以毫不夸张地说，大学本科高等数学中的大多数微积分问题，都能用符号计算解决，手工笔算演绎的烦劳都可以由计算机完成。
.3.1极限和导数的符号计算
limit,diff
【例2.3-1】两种重要极限
和
。
symstxk
s=sin（k*t）/（k*t）;
f=（1-1/x）^（k*x）;
Lsk=limit（s,0）%
Ls1=subs（Lsk,k,1）%
Lf=limit（f,x,inf）%
Lf1=vpa（subs（Lf,k,sym（'-1'））,48）%
Lsk=
1
Ls1=
1
Lf=
exp（-k）
Lf1=
2.7182818284590452353602874713526624977572470937
【例2.3-2】已知
，求
、
、
。
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]
【例2.3-3】求
的Jacobian矩阵
。
symsx1x2;
f=[x1*exp（x2）;x2;cos（x1）*sin（x2）];%
v=[x1;x2];%<3>
Jf=jacobian（f,v）%
Jf=
[exp（x2）,x1*exp（x2）]
[0,1]
[-sin（x1）*sin（x2）,cos（x1）*cos（x2）]
●指令<3>把自变量x1和x2写成列向量v。
但若写成v=[x1,x2]，所得Jacobian矩阵完全一样。
.3.2序列/级数的符号求和
S=symsum（f,v,a,b）即
v为整数
【例2.3-7】
1）
symsnk
f1=1/（k*（k+1））;
s1=symsum（f1,k,1,n）
s1=
n/（n+1）
.3.3符号积分
int（f,v）不定积分
int（f,v,a,b）定积分
【例2.3-8】
clear
symsabx
f=x*log（x）
s=int（f,x）
f=
x*log（x）
s=
（x^2*（log（x）-1/2））/2
【例2.3-10】
。
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=
（14912*2^（1/4））/4641-（6072064*2^（1/2））/348075+（64*2^（3/4））/225+1610027357/6563700
VF2=
224.92153573331143159790710032805
.4微分方程的符号解法
.4.1符号解法和数值解法的互补作用
.4.2求微分方程符号解的一般指令
dsolve（‘eq1,eq2,…’,‘cond1,cond2,…’,‘v’）
.4.3微分方程符号解示例
【例2.4-1】求
的解。
clearall%<1>
S=dsolve（'Dx=y,Dy=-x'）%
disp（''）
disp（['微分方程的解',blanks
（2）,'x',blanks（22）,'y']）
disp（[S.x,S.y]）
S=
x:
[1x1sym]
y:
[1x1sym]
微分方程的解xy
[C1*cos（t）+C2*sin（t）,C2*cos（t）–C1*sin（t）]
〖说明〗
【例2.4-2】微分方程
clearall%<1>
y=dsolve（'（Dy）^2-x*Dy+y=0','x'）%<2>
y=
x^2/4
-C5^2+x*C5
【例2.4-3】
。
1）
y=dsolve（'x*D2y-3*Dy=x^2','y
（1）=0,y（5）=0','x'）
y=
（31*x^4）/468-x^3/3+125/468
.5符号变换和符号卷积
.5.1Fourier变换及其反变换
（2.5-1）
（2.5-2）
Fw=fourier（ft,t,w）
ft=ifourier（Fw,w,t）
【例2.5-1】求单位阶跃函数的Fourier变换。
1）
symstw
ut=heaviside（t）;
UT=fourier（ut）
UT=
pi*dirac（w）-i/w
2）
Ut=ifourier（UT,w,t）%
SUt=simplify（Ut）%<5>
Ut=
（pi+pi*（2*heaviside（t）-1））/（2*pi）
SUt=
heaviside（t）
【例2.5-2】矩形脉冲
的Fourier变换。
1）
symsAtwtao
yt=A*（heaviside（t+tao/2）-heaviside（t-tao/2））;%
Yw=fourier（yt,t,w）%
Yws=simplify（Yw）%<4>
Yw=
A*（（cos（（tao*w）/2）*i+sin（（tao*w）/2））/w-（cos（（tao*w）/2）*i-sin（（tao*w）/2））/w）
Yws=
（2*A*sin（（tao*w）/2））/w
2）
Yt=ifourier（Yws,w,t）%
Yts=simplify（Yt）%<6>
Yt=
-（A*（pi*heaviside（t-tao/2）-pi*heaviside（t+tao/2）））/pi
Yts=
-A*（heaviside（t-tao/2）-heaviside（t+tao/2））
3）
T=3;
tn=-3:
0.1:
T;%
ytn=subs（yt,{A,tao,t},{1,T,tn}）;%
kk=find（tn==-t3/2|tn==t3/2）;%<11>
plot（tn（kk）,ytn（kk）,'.r','MarkerSize',30）%
ytn（kk）=NaN;%<13>
holdon
plot（tn,ytn,'-r','LineWidth',3）
holdoff
gridon
axis（[-T,T,-0.5,1.5]）
yt13=
heaviside（t+3/2）-heaviside（t-3/2）
图2.5-2矩形波
4）频域曲线绘制
Ywn=subs（Yws,{A,tao},{1,T}）;
subplot（2,1,1）,ezplot（Ywn）,gridon
subplot（2,1,2）,ezplot（abs（Ywn））,gridon
图2.5-3矩形脉冲的频率曲线和幅度频谱
.5.2Laplace变换及其反变换
（2.5-3）
（2.5-4）
Fs=laplace（ft,t,s）
ft=ilaplace（Fs,s,t）
【例2.5-4】求
的Laplace变换。
symsbpositive%<7>
f5=dirac（t-b）;%<11>
F5=laplace（f5,t,s）%
ft_F5=ilaplace（F5,s,t）%
F5=
exp（-b*s）
ft_F5=
dirac（b-t）
.5.3Z变换及其反变换
（2.5-5）
（2.5-6）
FZ=ztrans（fn,n,z）
fn=iztrans（FZ,z,n）
【例2.5-5】一组Z变换、反变换算例。
symsnzclear%<11>
f1=1;
F1=ztrans（f1,n,z）;
pretty（F1）
inv_F1=iztrans（F1,z,n）
z
-----
z-1
inv_F1=
1
.5.4符号卷积
【例2.5-6】已知系统冲激响应
，求
输入下的输出响应。
1）
clearall
symsTttaos
ut=exp（-t）;%
ht=exp（-t/T）/T;%
uh_tao=subs（ut,t,tao）*subs（ht,t,t-tao）;%
yt1=int（uh_tao,tao,0,t）%
yt1=
-（exp（-t）-exp（-t/T））/（T-1）
2）
yt21=ilaplace（laplace（ut,t,s）*laplace（ht,t,s）,s,t）
yt22=collect（yt21,'（T-1）'）%为与此前结果比较而施行的指令
yt21=
exp（-t/T）/（T-1）-exp（-t）/（T-1）
yt22=
（exp（-t/T）-exp（-t））/（T-1）
.6符号矩阵分析和代数方程解
.6.1符号矩阵分析
表2.6-1符号矩阵分析和解算指令汇总表
指令
含义
参考节次
colspace（A）
det（A）
diag（A）
[V,D]=eig（A）
expm（A）
inv（A）
[V,J]=jordan（A）
null（A）
charpoly（A）
rank（A）
rref（A）
s=svd（A）
[U,S,V]=svd（vpa（A））
tril（A）
triu（A）
〖说明〗
【例2.6-1】求矩阵
的行列式、逆和特征根。
1）
A=sym（'A%d%d',[2,2]）%
nA=ndims（A）%
SA=size（A）%
DA=det（A）
A=
[A11,A12]
[A21,A22]
nA=
2
SA=
22
DA=
A11*A22-A12*A21
2）
IA=inv（A）%
SIA=subexpr（IA,'ID'）%
IA=
[A22/（A11*A22-A12*A21）,-A12/（A11*A22-A12*A21）]
[-A21/（A11*A22-A12*A21）,A11/（A11*A22-
升级会员 