matlab ch2符号计算.docx
《matlab ch2符号计算.docx》由会员分享，可在线阅读，更多相关《matlab ch2符号计算.docx（24页珍藏版）》请在冰豆网上搜索。
matlabch2符号计算
第2章符号计算
符号计算：
解算数学表达式、方程不是在离散化的数值点上进行，而是凭借一系列恒等式，数学定理，通过推理和演绎，获得解析结果。
符号计算建立在数值完全准确表达和推演严格解析的基础之上，所得结果完全准确。
特点：
一，相对于MATLAB的数值计算“引擎”和“函数库”而言，符号计算的“引擎”和“函数库”是独立的。
二，在相当一些场合，符号计算解算问题的指令和过程，显得比数值计算更自然、更简明。
三，大多数理工科的本科学生在学过高等数学和其他专业基础课以后，比较习惯符号计算的解题理念和模式。
运算引擎MuPAD
MuPAD是极佳数学及符号数值运算绘图软件，同时也作为MATLAB7.8的符号计算工具箱，是一具有人工智能的数学软件,非常适合科学家及工程师使用.更适合每一个人使用,使用的方法非常简单,只要输入方程式就立刻得到答案,可以求Symbolic符号解,多项式之根,求非线性方程式之根,矩阵及向量VectorandMatrices运算,代数Algebra运算,求积分之值,求微分之值Calculus微积分等。
方程式可以处理复数计算.完美的绘图功能,图型输入,输出,轻松无比的绘图,可以输入多个2-D函数或极坐标函数或3-D函数,选择所要绘图参数,就可以完成图形，以及图形的动画制作也是非常方便。
数值计算结果并不是MATLAB命令行窗口所得的类似代码形式，而是规范数学格式。
并拥有一内建的程序语言，帮助文档以及文本操作，文本操作在一定程度上可以取代word.是一个超级的工程数学计算器.
MathWorks自从2008年10开始，在Matlab的新版本（Matlab2008a，即7.6之后）中使用MuPAD内核替换原来的Maple符号计算内核！
.1符号对象和符号表达式
MATLAB依靠基本符号对象（包括数字、参数、变量）、运算符及一些预定义函数来构造和衍生符号表达式和符号方程。
.1.1符号对象的创建和衍生
10一生成符号对象的基本规则
●任何基本符号对象都必须借助专门的符号函数指令sym或syms定义。
●任何包含符号对象的表达式或方程，将继承符号对象的属性。
10二符号数字
符号（类）数字的定义：
sym（'Num'）创建一个符号数字Num
sc=sym（'Num'）创建一个符号常数sc，该常数值准确等于Num
说明：
Num代表一个具体的数字
Num必须处于（英文状态下的）单引号内，构成字符串（关于字符串参见附录A.1）。
【例2.1-1】符号（类）数字与数值（类）数字之间的差异。
a=pi+sqrt（5）%创建方式
sa=sym（'pi+sqrt（5）'）
Ca=class（a）%类别判断
Csa=class（sa）
vpa（sa-a）
digits
digits（8）
10三符号参数
表达式e-axsinbx中的a,b称为参数。
定义格式：
symsPara定义符号参数Para
Para=sym（'Para'）
symsParaFlag定义具有Flag指定属性的符号参数Para
Para=sym（'Para','Flag'）
symsPara1Para2ParaN定义Para1Para2ParaN为符号参数
symsPara1Para2ParaNFlag定义Para1Para2ParaN为具有Flag指定属性的符号参数
●符号参数名不要用处于“字母表中小写字母x及其两侧的英文字母”开头。
●Flag表示参数属性，可具体取以下词条：
positive表示那些符号参数取正实数；
real表示那些符号参数限定为实时；
unreal/clear表示那些符号参数为不限定的复数。
10四符号变量
e-axsinbx中的x称为变量，符号变量的定义同符号参数。
确定自由符号变量的规则：
●在专门指定变量名的符号运算中，解题一定围绕指定变量名进行。
●自动识别符号变量时，字母的优先次序为x,y,w,z,v等。
自动识别表达式中自由、独立的符号变量的指令：
findsym（EXPR）确认表达式EXPR中所有自由符号变量
findsym（EXPR,N）确认表达式EXPR中距离x最近的N个自由符号变量
【例2.1-2】用符号计算研究方程
的解。
（1）不指定变量情况
symsuvwz%定义符号参数/变量
Eq=u*z^2+v*z+w;
result_1=solve（Eq）%
findsym（Eq,1）
（2）指定变量情况
result_2=solve（Eq,z）
【例2.1-3】对独立自由符号变量的自动辨认。
（1）
symsabxXY%定义符号参数/变量
k=sym（'3'）;%符号常数
z=sym（'c*sqrt（delta）+y*sin（theta1）'）;%直接定义符号表达式
EXPR=a*z*X+（b*x^2+k）*Y;%构成衍生符号表达式
（2）
findsym（EXPR）
（3）
findsym（EXPR,1）
（4）
findsym（EXPR,2）,findsym（EXPR,9）
【例2.1-4】findsym确定自由变量是对整个矩阵进行的。
symsabtuvxy
A=[a+b*x,sin（t）+u;x*exp（-t）,log（y）+v]
findsym（A,5）
A=
[a+b*x,u+sin（t）]
[x/exp（t）,v+log（y）]
ans=
x,y,v,u,t
.1.2符号计算中的算符
●与数值计算中的算符在形状、名称和使用方法上几乎完全相同。
●仅注意：
在符号对象的关系运算符中，只有算符“==”，“~=”
比较结果为“真”时，用1表示;否则用0表示。
.1.3符号计算中的函数指令
表2.1-1MATLAB中可调用的符号计算函数指令
类别
情况描述
与数值计算对应关系
基本函数
三角函数、双曲函数及反函数；除atan2外
名称和使用方法相同
指数、对数函数（如exp，expm）
名称和使用方法相同
复数函数（注意：
没有幅角函数angle）
名称和使用方法相同
矩阵分解函数（如eig，svd）
名称和使用方法相同
方程求解函数solve
不同
微积分函数（如diff，int）
不完全相同
积分变换和反变换函数（如laplace，ilaplace）
只有离散Fourier变换
绘图函数（如ezplot，ezsurf）
数值绘图指令更丰富
经典特殊函数
如误差函数erf、贝塞尔函数besselj、第一类完全椭圆积分EllipticK等；通过mfunlist可以看到所有经典函数名
mfunfunctionsperformnumerical,notsymbolic,calculations.Theinputparametersshouldbescalars,vectors,ormatricesoftypedouble,orcomplexdoubles,notsymbolicvariables.
部分
Maple库函数
Maple库函数在符号计算的扩展目录上；可通过mhelpindex看到各子函数库的名称；函数的数量很大；使用库函数，需要具备Maple语言知识
注意：
使用函数注意数据类型。
就数字而言，有双精度和符号类数字之分。
.1.4符号对象的识别
为了函数指令与数据对象的适配，MATLAB提供了用于识别数据对象属性的指令：
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;%产生4个数值变量
Mn=[a,b;c,d]%利用已赋值变量构成数值矩阵
Mc='[a,b;c,d]'%字符串中的a,b,c,d与前面输入的数值变量无关
Ms=sym（Mc）%Ms是一个符号矩阵，它与前面各变量无关
（2）
SizeMn=size（Mn）
SizeMc=size（Mc）
SizeMs=size（Ms）
（3）
CMn=class（Mn）
CMc=class（Mc）
CMs=class（Ms）
（4）
isa（Mn,'double'）
isa（Mc,'char'）
isa（Ms,'sym'）
（5）
whosMnMcMs
.2符号数字及表达式的操作
.2.1数值数字与符号数字之间的转换
10一数值数字向符号数字的转换
在符号运算中，“数值类数字”会自动地转换为符号数字。
亦可借助sym函数：
sym（Num,'r'）或sym（Num）数值类数字Num的广义有理表达
sym（Num,'d'）数值类数字Num的“十进制浮点”近似表达
sym（Num,'e'）数值类数字Num的带eps误差的理性近似表达
sym（Num,'f'）数值类数字Num的“十六进制浮点”近似表达
sym（'Num'）和sym（Num）区别问题
IEEE-754标准的浮点数格式
10二符号数字向双精度数字转换
double（Num_sym）把符号数字Num_sym转换为双精度数字
.2.2符号数字的任意精度计算
digits显示当前环境下符号数字“十进制浮点”表示的有效数字位数
digits（n）设定符号数字“十进制浮点”表示的有效数字位数（默认32位）
xs=vpa（x）据表达式x得到digits指定精度下的符号数字xs
xs=vpa（x,n）据表达式x得到n位有效数字的符号数字xs
【例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.3符号表达式的基本操作
collect（合并同类项）
factor（进行因式分解）
numden（提取公因式）等
最常用：
simple（EXPR）把EXPR转换成最简形式
【例2.2-2】简化
。
symsx
f=（1/x^3+6/x^2+12/x+8）^（1/3）;
g1=simple（f）
g2=simple（g1）
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.4表达式中的置换操作
10一子表达式置换操作
[RS,ssub]=subexpr（S,ssub）运用符号变量ssub置换子表达式，并重写S为RS
【例2.2-3】对符号矩阵
进行特征向量分解。
clearall
symsabcdW
[V,D]=eig（[ab;cd]）%V：
特征向量阵D：
特征值阵
[RVD,W]=subexpr（[V;D],W）%对矩阵元素中的公共子表达式进行置换表达
10二通用置换指令
RES=subs（ES,old,new）用new置换ES中的old后产生RES
RES=subs（ES,new）用new置换ES中的自由变量后产生RES
【例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）
（3）符号常数置换
f2=subs（f,{a,x},{2,sym（'pi/3'）}）%<3>a被双精度数字置换，x被符号数字置换
class（f2）
（4）双精度数值置换
f3=subs（f,{a,x},{2,pi/3}）%<4>
class（f3）
（5）数值数组置换之一
f4=subs（subs（f,a,2）,x,0:
pi/6:
pi）%<5>
class（f4）
（6）数值数组置换之二
f5=subs（f,{a,x},{0:
6,0:
pi/6:
pi}）%<6>
class（f5）
.3符号微积分
.3.1极限和导数的符号计算
大学本科高等数学中的大多数微积分问题，都能用符号计算解决，手工笔算演绎的烦劳都可以由计算机完成。
limit（f,v,a）求极限
【例2.3-1】试求
。
symsxk
Lim_f=limit（（1-1/x）^（k*x）,x,inf）
diff（f,v,n）求
（n缺省时，默认n=1）
【例2.3-2】求
求
，
，
。
symsatx;f=[a,t^3;t*cos（x）,log（x）];
df=diff（f）f对x的导数
dfdt2=diff（f,t,2）f对x的二阶导数
dfdxdt=diff（diff（f,x）,t）二阶混合导数
【例2.3-3】求
的Jacobian（雅可比）矩阵
。
symsx1x2;f=[x1*exp（x2）;x2;cos（x1）*sin（x2）];
v=[x1x2];fjac=jacobian（f,v）
.3.2序列/级数的符号求和
symsum（f,v,a,b）求f在变量v取遍[a,b]中所有整数时的和。
a,b缺省时默认求和区间[0,v-1]。
【例2.3-8】求
，
。
symskt;f1=[tk^3];f2=[1/（2*k-1）^2,（-1）^k/k];
s1=simple（symsum（f1））%f1的自变量被确认为t
s2=simple（symsum（f2,1,inf））%f1的自变量被确认为k
symsxy;f1=[yx^3];f2=[1/（2*y-1）^2,（-1）^y/y];
s1=simple（symsum（f1））%f1的自变量被确认为x
s2=simple（symsum（f2,1,inf））%f1的自变量被确认为y
.3.3符号积分
int（f,v）求f对变量v的不定积分
int（f,v,a,b）求f对变量v的定积分
【例2.3-9】求
。
clear
symsx
f=sqrt（（1+x）/x）/x
s=int（f,x）
s=simple（simple（s））
f=
（（x+1）/x）^（1/2）/x
s=
-2*（1/x+1）^（1/2）-2*atan（（1/x+1）^（1/2）*i）*i
s=
-2*（1/x+1）^（1/2）-2*atan（（1/x+1）^（1/2）*i）*i
【例2.3-10】求
。
symsabx;f=[a*x,b*x^2;1/x,sin（x）];
disp（'Theintegraloffis'）;pretty（int（f））
【例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位数字表示
【例2.3-12】求阿基米德（Archimedes）螺线
在
到
间的曲线长度函数，并求出
时的曲线长度。
（1）
symsartheta1phi1positive
x=r*cos（theta1）;x=subs（x,r,a*theta1）;
y=r*sin（theta1）;y=subs（y,r,a*theta1）;
dLdth=sqrt（diff（x,theta1）^2+diff（y,theta1）^2）;
L=simple（int（dLdth,theta1,0,phi1））
（2）
L_2pi=subs（L,[a,phi1],sym（'[1,2*pi]'））
L_2pi_vpa=vpa（L_2pi）
（3）
L1=subs（L,a,sym（'1'））;
ezplot（L1*cos（phi1）,L1*sin（phi1）,[0,2*pi]）
gridon
holdon
x1=subs（x,a,sym（'1'））;
y1=subs（y,a,sym（'1'））;
h1=ezplot（x1,y1,[0,2*pi]）;
set（h1,'Color','r','LineWidth',5）
title（''）
legend（'螺线长度-幅角曲线','阿基米德螺线'）
holdoff
1．说出以下三条指令产生的结果各属于哪种数据类型，是“双精度”对象，还是“符号”对象？
3/7+0.1,sym（3/7+0.1）,vpa（sym（3/7+0.1））
〖答案〗
c1=
0.5286
c2=
37/70
c3=
.52857142857142857142857142857143
2．在不加专门指定的情况下，以下符号表达式中的哪一个变量被认为是独立自由变量。
sym（'sin（w*t）'）,sym（'a*exp（-X）'）,sym（'z*exp（j*theta）'）
〖答案〗
ans=
w
ans=
a
ans=
z
3．求符号矩阵
的行列式值和逆，所得结果应采用“子表达式置换”简洁化。
〖答案〗
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
4．对函数
，当
为正实数时，求
。
（实际上，这就是根据定义求Z变换问题。
）
〖答案〗
Z1=
-z/（-z+a）
5．对于
，求
。
（提示：
理论结果为
）
〖答案〗
s_ss=
log（x）
6．
（1）通过符号计算求
的导数
。
（2）然后根据此结果，求
和
。
〖答案〗
d=
abs（1,sin（t））*cos（t）
d0_=
-1
dpi_2=
0
7．求出
的具有64位有效数字的积分值。
〖答案〗
matlab2008a
si=
1.087849499412904913166671875948174520895458535212845987519414167
8．计算二重积分
。
〖答案〗
r=
1006/105
9．在
区间，画出
曲线，并计算
。
〖答案〗
y=
sinint（x）
y5=
1.6541404143792439835039224868515
10．求
的一般积分表达式，并计算
的32位有效数字表达。
〖答案〗
yn=
1/2*pi^（1/2）*gamma（1/2+1/2*n）/gamma（1+1/2*n）
11．有序列
，
，（在此
，
），求这两个序列的卷积
。
〖答案〗
y1=
（-（b/a）^（k+1）+1）*a^k*a/（-b+a）
或
y2=
（-b*b^k+a^k*a）/（-b+a）
12．设系统的冲激响应为
，求该系统在输入
，
作用下的输出。
（提示：
运用卷积进行计算）
〖答案〗
hut=
-3/10/exp（t）^3+3/10*cos（t）+1/10*sin（t）
13．求
的Fourier变换。
〖答案〗
F=
2*A*a/（a^2+w^2）
14．求
的Fourier变换，并画出
时的幅频谱。
〖答案〗
Fws=
4*A/tao/w^2*sin（1/2*tao*w）^2
Fw2=
4/w^2*sin（w）^2
15．求
的Laplace反变换。
〖答案〗
f=
1/3*exp（-t）*（-2*cos（3^（1/2）*t）+3^（1/2）*sin（3^（1/2）*t）+2）
16．利用符号运算证明Laplace变换的时域求导性质：
。
〖答案〗
Ldy=
s*laplace（f（t）,t,s）-f（0）
17．求
的Z变换表达式。
〖答案〗
F_z=
z*exp（-lambda*T）/（z-exp（-lambda*T））^2
18．求方程
的解。
〖答案〗
x=
-1/2*（1/2*5^（1/2）+1/2*i*3^（1/2））^3+1/4*5^（1/2）+1/4*i*3^（1/2）
-1/2*（1/2*5^（1/2）-1/2*i*3^（1/2））^3+1/4*5^（1/2）-1/4*i*3^（1/2）
-1/2*（-1/2*5^（1/2）+1/2*i*3^（1/2））^3-1/4*5^（1/2）+1/4
升级会员 