matlab ch2例题doc文档格式.docx
《matlab ch2例题doc文档格式.docx》由会员分享，可在线阅读，更多相关《matlab ch2例题doc文档格式.docx（25页珍藏版）》请在冰点文库上搜索。
x,y,theta
【例2.1-4】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]
【例2.1-5】数据对象及其识别指令的使用。
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=
c,d]
Ms=
[a,b]
[c,d]
SizeMn=size（Mn）
SizeMc=size（Mc）
SizeMs=size（Ms）
SizeMn=
22
SizeMc=
19
SizeMs=
22
CMn=class（Mn）
CMc=class（Mc）
CMs=class（Ms）
CMn=
CMc=
char
CMs=
sym
isa（Mn,'
double'
）
isa（Mc,'
char'
isa（Ms,'
sym'
）
1
1
（5）
whosMnMcMs
NameSizeBytesClassAttributes
Mc1x918char
Mn2x232double
Ms2x2312sym
【例2.2-1】digits,vpa指令的使用。
digits
p0=sym（'
（1+sqrt（5））/2'
）
pr=sym（（1+sqrt（5））/2）
%
pd=sym（（1+sqrt（5））/2,'
d'
e32r=vpa（abs（p0-pr））
e16=vpa（abs（p0-pd）,16）
e32d=vpa（abs（p0-pd））
Digits=32
p0=
（1+sqrt（5））/2
pr=
7286977268806824*2^（-52）
pd=
1.6180339887498949025257388711907
e32r=
.543211520368251e-16
e16=
0.
e32d=
.543211520368251e-16
【例2.2-2】简化
。
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
【例2.2-3】对符号矩阵
进行特征向量分解。
clearall
symsabcdW
[V,D]=eig（[ab;
cd]）
[RVD,W]=subexpr（[V;
D],W）
V=
[-（1/2*d-1/2*a-1/2*（d^2-2*a*d+a^2+4*b*c）^（1/2））/c,-（1/2*d-1/2*a+1/2*（d^2-2*a*d+a^2+4*b*c）^（1/2））/c]
[1,1]
D=
[1/2*d+1/2*a+1/2*（d^2-2*a*d+a^2+4*b*c）^（1/2）,0]
[0,1/2*d+1/2*a-1/2*（d^2-2*a*d+a^2+4*b*c）^（1/2）]
RVD=
[-（1/2*d-1/2*a-1/2*W）/c,-（1/2*d-1/2*a+1/2*W）/c]
[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）
数学背景：
设A为n阶方阵,若存在数λ和非零的n维列向量x,使得Ax=λx或（A-λI）x=0，则称数λ为矩阵A的特征值,称x为矩阵A对应于特征值λ的特征向量.
【例2.2-4】用简单算例演示subs的置换规则。
（1）产生符号函数
symsax;
f=a*sin（x）+5
f=
a*sin（x）+5
（2）符号表达式置换
f1=subs（f,'
sin（x）'
sym（'
y'
））%<
2>
class（f1）
f1=
a*y+5
（3）ａ被双精度数字置换，ｘ被符号数字置换
f2=subs（f,{a,x},{2,sym（'
pi/3'
）}）%<
3>
class（f2）
f2=
3^（1/2）+5
sym
（4）ａ，ｘ均被双精度数值取代
f3=subs（f,{a,x},{2,pi/3}）%<
4>
class（f3）
f3=
6.7321
double
（5）ｘ被数值数组取代
f4=subs（subs（f,a,2）,x,0:
pi/6:
pi）%<
5>
class（f4）
f4=
5.00006.00006.73217.00006.73216.00005.0000
（6）ａ，ｘ被同样大小的数值数组取代，这时执行的是数组乘
f5=subs（f,{a,x},{0:
6,0:
pi}）%<
6>
class（f5）
f5=
5.00005.50006.73218.00008.46417.50005.0000
表达式能和符号变量同时包含在old中被置换，如
subs（f,{a,sin（x）},{2,sym（'
）}）
2*y+5
【例2.3-1】试求
symsxk
Lim_f=limit（（1-1/x）^（k*x）,x,inf）
Lim_f=
exp（-k）
【例2.3-2】求
求
，
symsatx;
f=[a,t^3;
t*cos（x）,log（x）];
df=diff（f）%求矩阵f对x的导数
dfdt2=diff（f,t,2）%求矩阵f对t的二阶导数
dfdxdt=diff（diff（f,x）,t）%求二阶混合导数
df=
[0,0]
[-t*sin（x）,1/x]
dfdt2=
[0,6*t]
dfdxdt=
[-sin（x）,0]
结论：
求导运算是对矩阵元素逐个进行的
【例2.3-3】求
的Jacobian矩阵
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）]
【例2.3-4】
，求
需要分区间求导
（1）对于x≥0,据导数定义求导
symsdpositive
f_p=sin（x）;
%x≥0时，f（x）改写而成
df_p=limit（（subs（f_p,x,x+d）-f_p）/d,d,0）%求x>
0区间的导数<
df_p0=limit（（subs（f_p,x,d）-subs（f_p,x,0））/d,d,0）%x=0+的右导数<
df_p=
cos（x）
df_p0=
1
（2）对于x≤0,据导数定义求导
f_n=sin（-x）;
%x≤0时，f（x）改写而成
df_n=limit（（f_n-subs（f_n,x,x-d））/d,d,0）%求x<
0区间的导数<
7>
df_n0=limit（（subs（f_n,x,0）-subs（f_n,x,-d））/d,d,0）%x=0-的左导数<
8>
df_n=
-cos（x）
df_n0=
-1
（3）直接利用diff求导
f=sin（abs（x））;
dfdx=diff（f,x） %只能得到ｘ≠０时的导函数
dfdx=
cos（abs（x））*abs（1,x）％abs（1,x）表示abs（x）的一阶导数，只在ｘ≠０时有值，等价于sign（x）.用abs（1,sym（'
x'
））可观察
（4）图形观察
xn=-3/2*pi:
pi/50:
0;
xp=0:
3/2*pi;
xnp=[xn,xp（2:
end）];
holdon
plot（xnp,subs（f,x,xnp）,'
k'
'
LineWidth'
2.5）%<
13>
plot（xn,subs（df_n,x,xn）,'
--r'
2.5）
plot（xp,subs（df_p,x,xp）,'
:
r'
2.5）
legend（char（f）,char（df_n）,char（df_p）,'
Location'
NorthEast'
）%<
16>
gridon
xlabel（'
holdoff
图2.3-1函数及其导函数
【例2.3-6】求
在
处展开的8阶Maclaurin级数。
r=taylor（x*exp（x）,9,x,0）%忽略9阶及9阶以上小量的展开
r=
x+x^2+1/2*x^3+1/6*x^4+1/24*x^5+1/120*x^6+1/720*x^7+1/5040*x^8
【例2.3-8】求
symskt;
f1=[tk^3];
f2=[1/（2*k-1）^2,（-1）^k/k];
findsym（f1,1）
findsym（f2）
s1=simple（symsum（f1））%f1的自变量被确认为t
s2=simple（symsum（f2,1,inf））%f2的自变量被确认为k
t
k
s1=
[1/2*t*（t-1）,k^3*t]
s2=
[1/8*pi^2,-log
（2）]
[说明]
●通式中的自变量只取整数值
●求和指令中的f可以是符号矩阵，此时求和操作将对矩阵中的“元素通式”逐个进行，但通式矩阵的自变量及其取值区间对各“元素通式”是相同的。
●由于findsym对f符号矩阵的自变量进行确认时，不是对该矩阵元素分别进行的，而是把矩阵f作为一个统一的对象处理的。
所以当f矩阵中含有一个以上非常数符号变量时要特别注意，求和对哪个符号变量进行。
【例2.3-9】求
f=sqrt（（1+x）/x）/x
s=int（f,x）%求不定积分
s=simple（simple（s））%简化计算结果
（（1+x）/x）^（1/2）/x
s=
（（1+x）/x）^（1/2）/x*（-2*（x^2+x）^（3/2）+2*（x^2+x）^（1/2）*x^2+log（1/2+x+（x^2+x）^（1/2））*x^2）/（（1+x）*x）^（1/2）
log（1/2+x+（（1+x）*x）^（1/2））-2*（（1+x）*x）^（1/2）/x
【例2.3-10】求
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）]
pretty指令把通常在一个物理行里显示的表达式扩展成用两个物理行显示的表达式。
目的是读起来更方便。
【例2.3-11】求积分
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）%积分结果用32位数字表示
F2=
1610027357/6563700+64/225*2^（3/4）-6072064/348075*2^（1/2）+14912/4641*2^（1/4）
VF2=
224.92153573331143159790710032805
对于内积分上下限为函数的多重积分，采用数值积分时编程较复杂。
【例2.3-12】求阿基米德（Archimedes）螺线
到
间的曲线长度函数，并求出
时的曲线长度。
symsarthetaphipositive%限定符号参数和变量取正值
x=r*cos（theta）;
x=subs（x,r,a*theta）;
y=r*sin（theta）;
y=subs（y,r,a*theta）;
dLdth=sqrt（diff（x,theta）^2+diff（y,theta）^2）;
L=simple（int（dLdth,theta,0,phi））
L=
1/2*a*（phi*（phi^2+1）^（1/2）-log（-phi+（phi^2+1）^（1/2）））
L_2pi=subs（L,[a,phi],sym（'
[1,2*pi]'
））%获得完全准确值
L_2pi_vpa=vpa（L_2pi）%计算32精度近似值
L_2pi=
pi*（1+4*pi^2）^（1/2）+1/2*asinh（2*pi）
L_2pi_vpa=
21.256294148209098800702512272566
L1=subs（L,a,sym（'
1'
））;
ezplot（L1*cos（phi）,L1*sin（phi）,[0,2*pi]）%ezplot的参变量使用格式
gridon%给坐标纸打方格线
holdon%保持上述图形窗，在其中继续绘图
x1=subs（x,a,sym（'
%使参数a数字化
y1=subs（y,a,sym（'
h1=ezplot（x1,y1,[0,2*pi]）;
%为改变ezplot绘制曲线色彩必须借助“图柄”
set（h1,'
Color'
5）%通过图柄改变曲线图形对象属性
legend（'
螺线长度-幅角曲线'
阿基米德螺线'
holdoff%在上述图形窗中，不再允许画任何图形
图2.3-2阿基米德螺线（粗红）和螺线长度函数（细蓝）
【例2.4-1】求
S=dsolve（'
Dx=y,Dy=-x'
）％从Ｓ的输出，只能看到２个域Ｓ.x和S.y
disp（[blanks（12）,'
blanks（21）,'
]）,disp（[S.x,S.y]）
S=
x:
[1x1sym]
y:
xy
[-C1*cos（t）+C2*sin（t）,C1*sin（t）+C2*cos（t）]
【例2.4-2】图示微分方程
的通解和奇解的关系。
y=dsolve（'
（Dy）^2-x*Dy+y=0'
）%注意书写规则，指定独立变量Ｘ
clf,holdon,
ezplot（y
（1）,[-6,6,-4,8],1）%画奇解
cc=get（gca,'
Children'
%取奇解曲线的图柄
set（cc,'
5）%把奇解画成粗红线
fork=-2:
0.5:
2;
ezplot（subs（y
（2）,'
C1'
k）,[-6,6,-4,8],1）;
end,%画通解，ezplot不能画多条曲线，必须循环解决
holdoff,title（'
\fontname{隶书}\fontsize{16}通解和奇解'
y=
1/4*x^2
-C1^2+x*C1
●dsolve可求通解和特解。
●若按题写成y=dsolve（'
y=x*Dy-（Dy）^2'
） 出错
图2.4-1通解和奇解曲线
【例2.5-1】求
的Fourier变换。
（1）求Fourier变换
symstw
ut=heaviside（t）;
%MATLAB7以前版本采用sym（'
Heaviside（t-a）'
UT=fourier（ut）
UT=
pi*dirac（w）-i/w
（2）求Fourier反变换
Ut=ifourier（UT,w,t）%结果应为原函数
Ut=
heaviside（t）
单位阶跃函数heaviside（t-a）=1,当t>
a;
heaviside（t-a）=0,当t<
a。
【例2.5-2】根据Fourier变换定义，用积分指令求方波脉冲
symsAtw
symstaopositive%限定是必须的
yt=heaviside（t+tao/2）-heaviside（t-tao/2）;
Yw=fourier（A*yt,t,w）
Yw=
2*A/w*sin（1/2*tao*w）
Yt=ifourier（Yw,w,t）
Yt=
A*（heaviside（t+1/2*tao）-heaviside（t-1/2*tao））
（3）绘制曲线
yt3=subs（yt,tao,3）
Yw3=subs（Yw,[A,tao],[1,3]）
subplot（2,1,1）
Ht=ezplot（yt3,[-3,3]）;
set（Ht,'
3）
subplot（2,1,2）,ezplot（Yw3）
yt3=
heaviside（t+3/2）-heaviside（t-3/2）
Yw3=
2/w*sin（3/2*w）
图2.5-1时域方波及其Fourier变换
【例2.5-3】求
的Fourier变换，在此
是参数，
是时间变量。
symstxw;
ft=exp（-（t-x））*heaviside（t-x）;
F1=simple（fourier（ft,t,w））%给出以w为频率变量的正确结果
F2=simple（fourier（ft））%误把x当作时间变量
F3=simple（fourier（ft,t））%误把x当作时间变量，又误把%t当作频率变量
F1=
1/（1+i*w）/exp（i*x*w）
i*exp（-i*t*w）/（i+w）
F3=
-exp（-t*（2+i*t））/（-1+i*t）
【例2.5-4】求
的Laplace变换。
本例演示：
符号参数的属性限定；
diric,heaviside的调用；
Laplace变换是对矩阵元素逐个实施的。
%positive
（1）在不对a,b限定时，运行结果可能不正确。
（2）对a，b进行限定后再运行可得正确结果
（3）再去除限定，运行仍然正确。
但当利用指令maplerestart重新启动Maple内核时，再在没限定情况下运行，仍将出错。
maplerestart
symsts;
symsab%positive
Dt=dirac（t-a）;
Ut=heaviside（t-b）;
Mt=[Dt,Ut;
exp（-a*t）*sin（b*t）,t^2*exp（-t）];
MS=laplace（Mt,t,s）
〖说明〗
●若不对符号参数a,b进行限定，不能确保heaviside（t-b）被正确变换。
●dirac（t-a）函数的数学含义是
在7.0版以前，它采用sym（'
Dirac（t-a）'
）实现。
比较而言，dirac（t-a）更便于使用。
【例2.5-6】求序列
的Z变换，并用反变换验算。
离散单位函数的构成；
ztrans,iztrans的用法。
（1）单位函数及性质验证
symsn
Delta=sym（'
charfcn[0]（n）'
%定义单位函数
<
D0=subs（Delta,n,0）;
%计算
D15=subs（Delta,n,15）;
%计算
[D0,D15]'
disp（[D0,D15]）
（2）求序列
的Z变换
symsz
fn=2*Delta+6*（1-（1/2）^n）%借助自定义的
升级会员 