matlab教程第6章.docx
- 文档编号:26891353
- 上传时间:2023-06-23
- 格式:DOCX
- 页数:24
- 大小:108.10KB
matlab教程第6章.docx
《matlab教程第6章.docx》由会员分享,可在线阅读,更多相关《matlab教程第6章.docx(24页珍藏版)》请在冰豆网上搜索。
matlab教程第6章
第6章符号计算
6.1符号对象和符号表达式
在MATLAB中,数值和数值变量用于数值的存贮和各种数值计算。
而符号常数、符号变量、符号函数、符号操作等则是用来形成符号表达式,严格按照代数、微积分等课程中的规则、公式进行运算,并尽可能给出解析结果。
6.1.1符号对象的生成和使用
f=sym(arg)把数字、字符串或表达式arg转换成为符号对象f
f=sym(argn,flagn)把数字、字符串或表达式argn转换成为flagn格式的符号对象f
flagn取:
‘d'——最接近的十进制浮点精确表示;
‘r'——最接近有理表示;
argv=sym('argv',flagv)按flagv指定的要求把字符串argv定义为符号对象argv
syms('argv1','argv2','argvk')把字符argv1,argv2,argvk定义为基本符号对象
symsargv1argv2argv3上述格式的简单形式
【例6.1-1】符号常数形成中的差异
a1=[1/3,pi/7,sqrt(5),pi+sqrt(5)]
a2=sym([1/3,pi/7,sqrt(5),pi+sqrt(5)])
a3=sym('[1/3,pi/7,sqrt(5),pi+sqrt(5)]')
a23=a2-a3
a1=
0.33330.44882.23615.3777
a2=
[1/3,pi/7,sqrt(5),6054707603575008*2^(-50)]
a3=
[1/3,pi/7,sqrt(5),pi+sqrt(5)]
a23=
[0,0,0,189209612611719/35184372088832-pi-5^(1/2)]
【例6.1-2】把字符表达式转换为符号变量
y=sym('2*sin(x)*cos(x)')
y=simple(y)
y=
2*sin(x)*cos(x)
y=
sin(2*x)
【例6.1-3】用符号计算验证三角等式
。
symsfai1fai2
y=simple(sin(fai1)*cos(fai2)-cos(fai1)*sin(fai2))
y=
sin(fai1-fai2)
【例6.1-4】求矩阵
的行列式值、逆和特征根
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)]
【例6.1-5】验证积分
。
symsAttaow
yf=int(A*exp(-i*w*t),t,-tao/2,tao/2);Yf=simple(yf)
Yf=
2*A*sin(1/2*tao*w)/w
6.2符号计算中的算符和基本函数
(1)基本运算符
矩阵算符+-*\/^
元素对元素算符.*./.\.^
6.2.1识别对象类别的指令
【例6.1-6】数据对象及其识别指令的使用。
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]
(1)判断矩阵的大小
SizeMn=size(Mn),SizeMc=size(Mc),SizeMs=size(Ms)
SizeMn=
22
SizeMc=
19
SizeMs=
22
(2)用class获得每种矩阵的类别
CMn=class(Mn),CMc=class(Mc),CMs=class(Ms)
CMn=
double
CMc=
char
CMs=
sym
(3)用isa判断矩阵的类别
isa(Mn,'double'),isa(Mc,'char'),isa(Ms,'sym')
ans=
1
ans=
1
ans=
1
(4)利用whos观察内存变量的类别和其它属性
whosMnMcMs
NameSizeBytesClass
Mc1x918chararray
Mn2x232doublearray
Ms2x2408symobject
Grandtotalis21elementsusing458bytes
6.2.2符号表达式中自由变量的确定
findsym(EXPR)确认表达式EXPR中所有"自由"符号"变量"。
findsym(EXPR,N)从表达式EXPR中确认出靠x最近的N个独立变量。
【例6.1-7】对独立自由符号变量的自动辨认。
symsabxXY;
k=sym('3');
z=sym('c*sqrt(delta)+y*sin(theta)');
EXPR=a*z*X+(b*x^2+k)*Y;
findsym(EXPR)
ans=
X,Y,a,b,c,delta,theta,x,y
findsym(EXPR,1)
ans=
x
findsym(EXPR,2),findsym(EXPR,3)
ans=
x,y
ans=
x,y,theta
【例6.1-8】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
6.3符号对象的操作和转换
6.3.1符号表达式的简化操作simple(expr)
【例6.2-1】简化
symsx
f=(1/x^3+6/x^2+12/x+8)^(1/3);
g1=simple(f)
g2=simple(g1)
g1=
(2*x+1)/x
g2=
2+1/x
6.3.2置换操作
6.3.2.1子表达式置换操作
[rsssub]=subexpr(S,ssub)运用符号变量ssub置换子表达式,重写S为rs
【例6.2-2】把复杂表达式中所含的多个相同子表达式用一个符号代替,使表达简洁。
clearall
symsabcdW
[V,D]=eig([ab;cd]);
[RVD,W]=subexpr([V;D],W)
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)
6.3.2.2通用置换指令
RES=subs(ES,old,new)用new置换ES中的old后产生RES
RES=subs(ES,new)用new置换ES中的自由变量后产生RES
【例6.2-3】用简单算例演示subs的置换规则。
symsax;f=a*sin(x)+5;
(1)符号变量置换
f1=subs(f,'sin(x)',sym('y'))
f1=
a*y+5
(2)符号常数置换
f2=subs(f,{a,x},{2,sym(pi/3)})
f2=
3^(1/2)+5
(3)双精度数值置换
f3=subs(f,{a,x},{2,pi/3})
f3=
6.7321
(4)数值数组置换一(取a=2,x=0:
pi/6:
pi)
f4=subs(subs(f,a,2),x,0:
pi/6:
pi)
f4=
5.00006.00006.73217.00006.73216.00005.0000
(4)数值数组置换二(取a=0:
6,x=0:
pi/6:
pi)
f5=subs(f,{a,x},{0:
6,0:
pi/6:
pi})
f5=
5.00005.50006.73218.00008.46417.50005.0000
6.3.2.3符号数值精度控制和任意精度计算
【例6.2-4】指令使用演示。
p0=sym('(1+sqrt(5))/2');
p1=sym((1+sqrt(5))/2)
pd=double(p0)
e01=vpa(abs(p0-p1))
e0d=vpa(abs(p0-pd))
p1=
7286977268806824*2^(-52)
pd=
1.6180
e01=
.543211520368250e-16
e0d=
.543211520368250e-16
p2=vpa(p0)
e02=vpa(abs(p0-p2),64)
p2=
1.6180339887498948482045868343657
e02=
.61882279690820194237137864551377e-31
digits
Digits=32
6.3.3符号对象与其它数据对象间的转换
6.4符号微积分
6.4.1符号序列的求和
s=symsum(f,v,a,b)求通式f在指定变量v取遍[a,b]中所有整数时的和
【例6.3-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)]
6.4.2符号微分和
矩阵
【例6.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]
【例6.3.2-2】求
的
矩阵。
symsx1x2;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)]
6.4.3符号积分
intf=int(f,v)给出f对指定变量v的不定积分
intf=int(f,v,a,b)给出f对指定变量v的定积分
【例6.3-4】求
。
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)]
【例6.3-5】求
。
F1=int('1/log(t)','t',0,'x')
F1=
-Ei(1,-log(x))
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
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=
1610027357/6563700-6072064/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2^(3/4)
VF2=
224.92153573331143159790710032805
6.4.4符号卷积
【例5.3-7】本例演示卷积的时域积分法:
已知系统冲激响应
,求
输入下的输出响应。
symsTtpositive
symstao
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)
【例6.3-8】本例演示通过变换和反变换求取卷积。
系统冲激响应、输入同上例,求输出。
symss
yt=ilaplace(laplace(ut,t,s)*laplace(ht,t,s),s,t);
yt=simple(yt)
yt=
-(exp(-t)-exp(-t/T))/(T-1)
【例6.3-9】求函数
和
的卷积。
symstao;
t=sym('t','positive');
ut=sym('Heaviside(t)-Heaviside(t-1)');ht=t*exp(-t);
yt53=int(subs(ut,t,tao)*subs(ht,t,t-tao),tao,0,t);
yt53=collect(yt53,'Heaviside(t-1)')
yt53=
(-1+exp(1-t)*t)*Heaviside(t-1)+1+(-t-1)*exp(-t)
6.4.5符号积分变换
6.4.5.1Fourier变换及其反变换
【例6.4-1】求
的Fourier变换。
symstw
ut=sym('Heaviside(t)');
UT=fourier(ut)
UTC=maple('convert',UT,'piecewise','w')
UTS=simple(UT)
UT=
pi*Dirac(w)-i/w
UTC=
PIECEWISE([undefined,w=0],[0,otherwise])
UTS=
-i/w
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)
【例6.4-2】用fourier指令求方波脉冲
的Fourier变换。
(1)
symsAtw
symstaopositive
yt=sym('Heaviside(t+tao/2)-Heaviside(t-tao/2)');
Yw=fourier(A*yt,t,w)
Ywc=maple('convert',Yw,'piecewise','w')
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))
【例6.4-3】求
的Fourier变换,在此
是参数,
是时间变量。
本例演示:
fourier的缺省调用格式的使用要十分谨慎;在被变换函数中包含多个符号变量的情况下,对被变换的自变量给予指明,可保证计算结果的正确。
symstxw;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)
6.4.5.2Laplace变换及其反变换
【例6.4-4】求
的Laplace变换。
symsts;symsabpositive
Dt=sym('Dirac(t-a)');
Ut=sym('Heaviside(t-b)');
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]
【例6.4-5】验证Laplace时移性质:
。
symsts;t0=sym('t0','positive');
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)
6.4.5.3Z变换及其反变换
【例6.4-6】求序列
的Z变换,并用反变换验算。
symsn
Delta=sym('charfcn[0](n)');
D0=subs(Delta,n,0);
D15=subs(Delta,n,15);
disp('[D0,D15]');disp([D0,D15])
[D0,D15]
[1,0]
symsz;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
4z+2
--------------
2
2z-3z+1
FZ_n=
2*charfcn[0](n)+6-6*(1/2)^n
6.5符号代数方程的求解
6.5.1线性方程组的符号解
【例6.5-1】求
线性方程组的解。
A=sym([11/21/2-1;11-11;1-1/4-11;-8-111]);
b=sym([0;10;0;1]);X1=A\b
X1=
[1]
[8]
[8]
[9]
6.5.2一般代数方程组的解
S=solve('eq1','eq2',…,'eqn','v1','v2',…,'vn')
【例6.5-2】求方程组
,
关于
的解。
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:
[2x1sym]
z:
[2x1sym]
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*w-v-(4*u*w*v+v^2-4*u*w)^(1/2))]
【例6.5-3】solve指令求
,
,
构成的“欠定”方程组解。
symsdnpq
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:
3equationsin4variables.
>InE:
\MAT53\toolbox\symbolic\solve.matline110
InE:
\MAT53\toolbox\symbolic\@sym\solve.matline49
ans=
d
ans=
8
ans=
4*d+4
ans=
3*d+6
【例6.5-4】求
的解。
clearall,symsx;s=solve('(x+2)^x=2','x')
s=
.69829942170241042826920133106081
6.6符号微分方程的求解
求微分方程符号解的一般指令
S=dsolve('a_1','a_2',…,'a_n')
6.6.1微分方程符号解示例
【例6.6-1】求
的解。
S=dsolve('Dx=y,Dy=-x');
disp([blanks(12),'x',blanks(21),'y']),disp([S.x,S.y])
xy
[cos(t)*C1+sin(t)*
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 教程