probans210a.docx
《probans210a.docx》由会员分享，可在线阅读，更多相关《probans210a.docx（25页珍藏版）》请在冰豆网上搜索。
probans210a
第2章符号运算
习题2及解答
11说出以下四条指令产生的结果各属于哪种数据类型，是“双精度”对象，还是“符号”符号对象？
3/7+0.1;sym（3/7+0.1）;sym（'3/7+0.1'）;vpa（sym（3/7+0.1））
〖目的〗
●不能从显示形式判断数据类型，而必须依靠class指令。
〖解答〗
c1=3/7+0.1
c2=sym（3/7+0.1）
c3=sym（'3/7+0.1'）
c4=vpa（sym（3/7+0.1））
Cs1=class（c1）
Cs2=class（c2）
Cs3=class（c3）
Cs4=class（c4）
c1=
0.5286
c2=
37/70
c3=
0.52857142857142857142857142857143
c4=
0.52857142857142857142857142857143
Cs1=
double
Cs2=
sym
Cs3=
sym
Cs4=
sym
12在不加专门指定的情况下，以下符号表达式中的哪一个变量被认为是自由符号变量.
sym（'sin（w*t）'）,sym（'a*exp（-X）'）,sym（'z*exp（j*th）'）
〖目的〗
●理解自由符号变量的确认规则。
〖解答〗
symvar（sym（'sin（w*t）'）,1）
ans=
w
symvar（sym（'a*exp（-X）'）,1）
ans=
a
symvar（sym（'z*exp（j*th）'）,1）
ans=
z
13求以下两个方程的解
（1）试写出求三阶方程
正实根的程序。
注意：
只要正实根，不要出现其他根。
（2）试求二阶方程
在
时的根。
〖目的〗
●体验变量限定假设的影响
〖解答〗
（1）求三阶方程
正实根
reset（symengine）%确保下面操作不受前面指令运作的影响
symsxpositive
solve（x^3-44.5）
ans=
（2^（2/3）*89^（1/3））/2
（2）求五阶方程
的实根
symsapositive%注意：
关于x的假设没有去除
solve（x^2-a*x+a^2）
Warning:
Explicitsolutioncouldnotbefound.
>Insolveat81
ans=
[emptysym]
symsxclear
symsapositive
solve（x^2-a*x+a^2）
ans=
a/2-（3^（1/2）*a*i）/2
a/2+（3^（1/2）*a*i）/2
14观察一个数（在此用@记述）在以下四条不同指令作用下的异同。
a=@,b=sym（@）,c=sym（@,'d'）,d=sym（'@'）
在此，@分别代表具体数值7/3,pi/3,pi*3^（1/3）；而异同通过vpa（abs（a-d））,vpa（abs（b-d））,vpa（abs（c-d））等来观察。
〖目的〗
●理解准确符号数值的创建法。
●高精度误差的观察。
〖解答〗
（1）x=7/3
x=7/3;a=x,b=sym（x）,c=sym（x,'d'）,d=sym（'7/3'）,
a=
2.3333
b=
7/3
c=
2.3333333333333334813630699500209
d=
7/3
v1=vpa（abs（a-d））,v2=vpa（abs（b-d））,v3=vpa（abs（c-d））
v1=
0.
v2=
0.
v3=
.1480297366166876e-15
（2）x=pi/3
x=pi/3;a=x,b=sym（x）,c=sym（x,'d'）,d=sym（'pi/3'）,
a=
1.0472
b=
pi/3
c=
.0471********
d=
pi/3
v1=vpa（abs（a-d））,v2=vpa（abs（b-d））,v3=vpa（abs（c-d））
v1=
0.
v2=
0.
v3=
.1148364282799222e-15
（3）x=pi*3^（1/3）
x=pi*3^（1/3）;a=x,b=sym（x）,c=sym（x,'d'）,d=sym（'pi*3^（1/3）'）
a=
4.5310
b=
5101408179057732*2^（-50）
c=
4.5309606547207899041040946030989
d=
pi*3^（1/3）
v1=vpa（abs（a-d））,v2=vpa（abs（b-d））,v3=vpa（abs（c-d））
v1=
.2660111416629094e-15
v2=
.2660111416629094e-15
v3=
.2660111416629094e-15
15求符号矩阵
的行列式值和逆，所得结果应采用“子表达式置换”简洁化。
〖目的〗
●理解subexpr指令。
〖解答〗
A=sym（'[a11a12a13;a21a22a23;a31a32a33]'）
DA=det（A）
IA=inv（A）;
[IAs,d]=subexpr（IA,d）
A=
[a11,a12,a13]
[a21,a22,a23]
[a31,a32,a33]
DA=
a11*a22*a33-a11*a23*a32-a21*a12*a33+a21*a13*a32+a31*a12*a23-a31*a13*a22
IAs=
[（a22*a33-a23*a32）/d,-（a12*a33-a13*a32）/d,-（-a12*a23+a13*a22）/d]
[-（a21*a33-a23*a31）/d,（a11*a33-a13*a31）/d,-（a11*a23-a13*a21）/d]
[-（-a21*a32+a22*a31）/d,-（a11*a32-a12*a31）/d,（a11*a22-a12*a21）/d]
d=
a11*a22*a33-a11*a23*a32-a21*a12*a33+a21*a13*a32+a31*a12*a23-a31*a13*a22
16对函数
，当
为正实数时，求
。
（实际上，这就是根据定义求Z变换问题）
〖目的〗
●symsum,ztrans,simple指令的应用。
〖解答〗
symskaz;
f=a^k*z^（-k）;
Z1=simple（symsum（f,k,0,inf））
disp（'用Z变换检验'）
Z2=simple（ztrans（a^k））%用Z变换检验
Z1=
-z/（-z+a）
用Z变换检验
Z2=
-z/（-z+a）
17对于
，求
。
（提示：
理论结果为
）
〖目的〗
●符号变量的限定性定义的作用。
〖解答〗
symsk;
x=sym（'x','positive'）;
f_k=2/（2*k+1）*（（x-1）/（x+1））^（2*k+1）;
s=simple（symsum（f_k,k,0,inf））%结果与理论值lnx相符！
s=
log（x）
18
（1）通过符号计算求
的导数
。
（2）然后根据此结果，求
和
。
〖目的〗
●diff,limit指令的应用。
●如何理解运行结果。
〖解答〗
法一：
symst
y=abs（sin（t））
d=diff（y）%求dy/dt
d0_=limit（d,t,0,'left'）%求dy/dt|t=0-
dpi_2=limit（d,t,pi/2）%求dy/dt|t=pi/2
y=
abs（sin（t））
d=
sign（sin（t））*cos（t）
d0_=
-1
dpi_2=
0
法二：
clear
symst
symsdpositive
f_p=sin（t）;%x≥0时，f（x）改写而成
df_p=limit（（subs（f_p,t,t+d）-f_p）/d,d,0）%求x>0区间的导数<4>
df_p0=subs（df_p,t,sym（‘pi/2’））%x=0+的右导数<5>
f_n=sin（-t）;%x≤0时，f（x）改写而成
df_n=limit（（f_n-subs（f_n,t,t-d））/d,d,0）%求x<0区间的导数<7>
df_n0=limit（（subs（f_n,t,0）-subs（f_n,t,-d））/d,d,0）%x=0-的左导数<8>
19求出
的具有64位有效数字的积分值。
〖目的〗
●符号积分的解析解和符号数值解。
●符号计算和数值计算的相互校验。
〖解答〗
（1）符号积分
y='exp（-abs（x））*abs（sin（x））'
si=vpa（int（y,-10*pi,1.7*pi）,64）
y=
exp（-abs（x））*abs（sin（x））
si=
1.087849499412904913166671875948174520895458535212845987519414166
（2）数值计算复验
xx=-10*pi:
pi/100:
1.7*pi;
sn=trapz（exp（-abs（xx））.*abs（sin（xx）））*pi/100
sn=
1.0877
110计算二重积分
。
〖目的〗
●变上限二重积分的符号计算法。
〖解答〗
symsxy
f=x^2+y^2;
r=int（int（f,y,1,x^2）,x,1,2）
r=
1006/105
111在
区间，画出
曲线，并计算
。
〖目的〗
●在符号计算中，经常遇到计算结果是特殊经典函数的情况。
●如何应用subs获得超过16位有效数字的符号数值结果。
●初步尝试ezplot指令的简便。
〖解答〗
（1）符号计算
symstx;
f=sin（t）/t;
y=int（f,t,0,x）%将得到一个特殊经典函数
y5=subs（y,x,sym（'4.5'））
ezplot（y,[0,2*pi]）
y=
sinint（x）
y5=
1.6541404143792439835039224868515
（2）数值计算复验
tt=0:
0.001:
4.5;
tt
（1）=eps;
yn=trapz（sin（tt）./tt）*0.001
yn=
1.6541
112在
的限制下，求
的一般积分表达式，并计算
的32位有效数字表达。
〖目的〗
●一般符号解与高精度符号数值解。
〖解答〗
symsx
symsnpositive
f=sin（x）^n;
yn=int（f,x,0,pi/2）
y3s=vpa（subs（yn,n,sym（'1/3'）））
y3d=vpa（subs（yn,n,1/3））
yn=
beta（1/2,n/2+1/2）/2
y3s=
1.2935547796148952674767575125656
y3d=
1.2935547796148951782413405453553
113有序列
，
，（在此
，
），求这两个序列的卷积
。
〖目的〗
●符号离散卷积直接法和变换法。
〖解答〗
（1）直接法
symsabkn
x=a^k;
h=b^k;
w=symsum（subs（h,k,n）*subs（x,k,k-n）,n,0,k）%据定义
y1=simple（w）
w=
piecewise（[a=b,b^k+b^k*k],[a<>b,（a*a^k-b*b^k）/（a-b）]）
y1=
piecewise（[a=b,b^k+b^k*k],[a<>b,（a*a^k-b*b^k）/（a-b）]）
（2）变换法（复验）
symsz
X=ztrans（a^k,k,z）;
H=ztrans（b^k,k,z）;
y2=iztrans（H*X,z,k）%通过Z变换及反变换
y2=
piecewise（[a<>0andb<>0,kroneckerDelta（k,0）+（a^2*（a^k/a-kroneckerDelta（k,0）/a））/（a-b）-（b^2*（b^k/b-kroneckerDelta（k,0）/b））/（a-b）]）
〖说明〗
●符号计算不同途径产生的结果在形式上有可能不同，而且往往无法依靠符号计算本身的指令是它们一致。
此时，必须通过手工解决。
114设系统的冲激响应为
，求该系统在输入
，
作用下的输出。
〖目的〗
●符号连续函数卷积的直接法和变换法。
●符号变量限定性定义的作用。
●laplace,ilaplace指令的应用。
〖解答〗
（1）直接法
symst
h=exp（-3*t）;u=cos（t）;
symstao;
h_tao=subs（h,t,tao）;
u_t_tao=subs（u,t,t-tao）;
hu_tao=h_tao*u_t_tao;
hut=simple（int（hu_tao,tao,0,t））%直接卷积
hut=
（3*cos（t））/10-3/（10*exp（3*t））+sin（t）/10
（2）变换法（复验）
symss;
HU=laplace（h,t,s）*laplace（u,t,s）;
huL=simple（ilaplace（HU,s,t））%拉氏变换及反变换
huL=
（3*cos（t））/10-3/（10*exp（3*t））+sin（t）/10
115求
的Fourier变换。
〖目的〗
●符号变量限定性定义的作用。
●fourier指令的应用。
〖解答〗
symsAtw
a=sym（'a','positive'）;
f=A*exp（-a*abs（t））;
y=fourier（f,t,w）
F=simple（y）
y=
（2*A*a）/（a^2+w^2）
F=
（2*A*a）/（a^2+w^2）
116求
的Fourier变换，并画出
时的幅频谱。
〖目的〗
●单位阶跃符号函数heaviside的应用。
●subs实现多变量置换。
●ezplot的使用。
〖解答〗
symstAw;
tao=sym（'tao','positive'）;
f=A*（（1+t/tao）*（heaviside（t+tao）-heaviside（t））+（1-t/tao）*（heaviside（t）-heaviside（t-tao）））;
Fw=fourier（f,t,w）;
Fws=simple（Fw）
Fw2=subs（Fws,[A,tao],[2,2]）
ezplot（abs（Fw2））
grid
Fws=
（4*A*sin（（tao*w）/2）^2）/（tao*w^2）
Fw2=
（4*sin（w）^2）/w^2
117求
的Laplace反变换。
〖解答〗
symsst
F=（s+3）/（s^3+3*s^2+6*s+4）;
f=simple（ilaplace（F,s,t））
f=
（3^（1/2）*sin（3^（1/2）*t）-2*cos（3^（1/2）*t）+2）/（3*exp（t））
118利用符号运算证明Laplace变换的时域求导性质：
。
〖目的〗
●符号计算用于定理证明。
〖解答〗
symsts;
y=sym（'f（t）'）;
df=diff（y,t）;
Ldy=laplace（df,t,s）
Ldy=
s*laplace（f（t）,t,s）-f（0）
119求
的Z变换表达式。
〖目的〗
●注意：
变换中，被变换变量的约定。
〖解答〗
symslambdakTz;
f_k=k*exp（-lambda*k*T）;
F_z=simple（ztrans（f_k,k,z））
F_z=
（z*exp（T*lambda））/（z*exp（T*lambda）-1）^2
120求方程
的解。
〖目的〗
●solve指令中，被解方程的正确书写，输出量的正确次序。
〖解答〗
eq1='x^2+y^2=1';
eq2='x*y=2';
[x,y]=solve（eq1,eq2,'x','y'）
x=
（1/2+（15^（1/2）*i）/2）^（1/2）/2-（1/2+（15^（1/2）*i）/2）^（3/2）/2
-（1/2+（15^（1/2）*i）/2）^（1/2）/2+（1/2+（15^（1/2）*i）/2）^（3/2）/2
（1/2-（15^（1/2）*i）/2）^（1/2）/2-（1/2-（15^（1/2）*i）/2）^（3/2）/2
-（1/2-（15^（1/2）*i）/2）^（1/2）/2+（1/2-（15^（1/2）*i）/2）^（3/2）/2
y=
（1/2+（15^（1/2）*i）/2）^（1/2）
-（（15^（1/2）*i）/2+1/2）^（1/2）
（1/2-（15^（1/2）*i）/2）^（1/2）
-（1/2-（15^（1/2）*i）/2）^（1/2）
121求图p2-1所示信号流图的系统传递函数，并对照胡寿松主编“自动控制原理”中的例2-21结果，进行局部性验证。
图p2-1
〖目的〗
●理解和掌握信号流图传递函数的“代数状态方程解法”。
●并设法用胡寿松主编的“自动控制原理”的例2-21进行局部性验证。
〖解答〗
（1）求传递函数
symsG1G2G3G4G5G6G7H1H2H3H4H5
A=[0000-H3-H4;
G10-H1000;
0G200-H2G6;
00G300G7;
000G400;
0G5000-H5];
b=[1;0;0;0;0;0];
c=[000010];
Y2U=c*（（eye（size（A））-A）\b）;%求传递函数
[NN,DD]=numden（Y2U）;%分离出分子、分母多项式
DD=sort（DD）;%分母多项式排序
disp（[blanks（5）,'传递函数Y2U为']）
pretty（NN/DD）
传递函数Y2U为
（G1G4（G2G3+G5G7+G3G5G6+G2G3H5））/
（H5+G2H1+G3G4H2+G1G5H4+G5G6H1+G2H1H5+G3G4H2H5+
G1G2G3G4H3+G1G4G5G7H3-G4G5G7H1H2+G1G3G4G5G6H3+
G1G2G3G4H3H5+G1G3G4G5H2H4+1）
（2）局部性验证
symsabcdefg
y2u=subs（Y2U,[G1,G2,G3,G4,G5,G6,G7,H1,H2,H3,H4,H5],[a,e,f,1,b,c,0,g,0,0,0,d]）;
[nn,dd]=numden（y2u）;
dd=sort（dd）;
disp（[blanks（5）,'局部性验证用的传递函数y2u']）
pretty（nn/dd）
局部性验证用的传递函数y2u
a（ef+bcf+def）
---------------------------
d+eg+bcg+deg+1
此结果与胡寿松主编的“自动控制原理”例2-21一致。
122采用代数状态方程法求图p2-2所示结构框图的传递函数
和
。
图p2-2
〖目的〗
●运用“代数状态方程解法”求输入和扰动同时存在的结构框图的传递函数。
〖解答〗
（1）理论演绎
对于结构框图写出状态方程
（p2-1）
此式第一个方程关于x的解可写为
（p2-2）
把此式代入式（p2-1）的第二个方程，加以整理后可得
据此可写出传递函数
（p2-3）
（p2-4）
（2）列出“元素级”状态方程
值得提醒：
在编写M码之前，最好先在草稿纸上，仔细“元素级”状态方程是避免出错的冲要措施。
对此，不要掉以轻心。
本例的“元素级”状态方程如下
（p2-5）
（3）编写相应的M码
symsG1G2G3H1H2
A=[000-G1-G1;
G20-G200;
00000;
0H1000;
0H2000];
b=[G1;0;0;0;0];
f=[0;0;G3;0;-H2];
c=[01000];
d=0;
g=-1;
R=c/（eye（size（A））-A）;%中间变量
Y2U=R*b+d;%计算传递函数Y/U
Y2W=R*f+g;%计算传递函数Y/W
[NU,DU]=numden（Y2U）;%分离出分子、分母多项式
DU=sort（DU）;%分母多项式排序
disp（[blanks（5）,'传递函数Y2U为']）
pretty（NU/DU）
[NW,DW]=numden（Y2W）;
NW=sort（NW）;
DW=sort（DW）;
disp（[blanks（5）,'传递函数Y2W为']）
pretty（NW/DW）
传递函数Y2U为
G1G2
-----------------------
G1G2H1+G1G2H2+1
传递函数Y2W为
G2G3+G1G2H1+1
------------------------
G1G2H1+G1G2H2+1
123求微分方程
的通解，并绘制任意常数为1时解的图形。
〖目的〗
●理解指令dsolve的正确
升级会员 