第3讲 matlab的符号运算.docx
《第3讲 matlab的符号运算.docx》由会员分享，可在线阅读，更多相关《第3讲 matlab的符号运算.docx（21页珍藏版）》请在冰点文库上搜索。
第3讲matlab的符号运算
第三讲MATLAB的符号运算
（注：
文中红色字体为命令执行的结果，在Command窗口中显示）
3-1符号对象的创建和使用
1．符号运算入门
符号运算的特点是，运算过程中允许存在非数值的符号变量。
先看如下示例：
函数
，用MATLAB求它的微积分，命令如下：
f=’sin（x）^2’;%定义符号函数f（x）
dfdx=diff（f）%求
的指令
intf=int（f）%求
的指令
显示的计算结果为：
dfdx=
2*sin（x）*cos（x）
intf=
-1/2sin（x）*cos（x）+1/2*x
所以，
，
。
此例中，首先定义符号函数f=’sin（x）^2’，然后由符号运算获得
的微分和积分。
2．定义符号变量
在使用符号变量之前，应先声明某些要用到的变量是“符号”变量。
声明符号变量的语句：
syms变量名列表
或：
sym（‘变量名’）
其中各个变量名应该用空格分隔，而不能用逗号分隔。
如创建符号变量x和a：
x=sym（‘x’）
a=sym（‘alpha’）
或用：
symsxa%定义符号变量x和a
这里，变量x和a的类型是符号对象，它们被定义后，即可参与符号运算。
3．定义符号表达式和符号方程
符号表达式和符号方程是两种不同的操作对象。
区别在于：
符号表达式不包含等号（=），而符号方程须带等号。
它们的创建方式相同。
如：
要考虑二次函数f=ax^2+bx+c，可以创建符号表达式，赋值给符号变量f。
f=sym（‘a*x^2+b*x+c’）
或：
f=‘a*x^2+b*x+c’
此例中，将符号表达式赋给符号变量f，但这不是必需的，引入符号变量是为了以后调用方便。
在这种情况下，没有创建对应于表达式中a、b、c、x项的变量，为了执行符号数学运算（如微分、积分等），必须显式地创建这些变量，可用下列命令创建：
symsabcx
如下例中创建了符号表达式和符号方程，分别赋给相应的符号对象。
symsxabc
f=’sin（x）^2’;%创建符号表达式sin（x）^2赋给变量f
eq=‘a*x^2+b*x+c=0’%创建的符号方程赋给变量eq
4．定义抽象函数和符号数学函数
若要创建抽象函数f（x），可使用：
f=sym（‘f（x）’）
则f就象f（x）一样参与运算。
在SymbolicMathToolbox中可利用符号表达式创建符号数学函数。
如下面的命令将产生符号表达式r、t、f：
symsxyz
r=sqrt（x^2+y^2+z^2）%或r=sym（‘sqrt（x^2+y^2+z^2）’）
t=atan（y/x）
f=sin（x*y）/（x*y）
3-2数值与符号的转换
Sym函数有四种选项将数值结果转换为符号表达式。
如对于变量rho，定义为：
，在MATLAB命令窗口中定义rho:
rho=（1+sqrt（5））/2
将这一数值结果转化为符号表达式：
（1）sym（rho,’f’）返回符号浮点表示形式，结果为：
sym（rho,’f’）
ans=
‘1.9e3779b97f4a8’*2^（0）
（2）sym（rho,’r’）返回符号有理数表示形式，这是sym的默认设置，当调用sym而没有第二个参数时，相当于sym使用了r选项。
结果为：
sym（rho,’r’）
ans=
7286977268806824*2^（-52）
例如：
sym（0.75）
ans=3/4
（3）sym（rho,’e’）返回符号有理数表示形式，同时根据eps给出rho的理论表达式和实际计算的差。
sym（rho,’e’）
ans=7286977268806824*2^（-52）
这里，rho的理论值和实际浮点值相同，对于1/3，可得：
sym（1/3,’e’）
ans=1/3-eps/12
（4）sym（rho,’d’）返回符号十进制小数表示形式。
有效位数由digits定义。
Digits的默认值是32位。
sym（rho,’d’）
ans=1.6180339887498949025257388711907
如果想要一个较短的结果，可由digits定义有效位数。
digits（5）
sym（rho,’d’）
ans=
1.6180
上述四种格式中，最后得到的都是符号对象变量。
3-3符号算术运算
1．定义符号矩阵
MATLAB中的数值矩阵不能直接参与符号运算，必须经过转换。
Sym函数的一个非常有用的功能就是将数值矩阵转换成符号矩阵。
例如，下面的代码将产生一个3×3的Hilbert矩阵。
A=hilb（3）
A=
1.00000.50000.3333
0.50000.33330.2500
0.33330.25000.2000
对A使用sym命令可获得一个具有无穷精度的3×3的符号形式的矩。
A=sym（A）
A=
[1,1/2,1/3]
[1/2,1/3,1/4]
[1/3,1/4,1/5]
当调用SymbolicMathToolbox时，默认地调用有理算术运算。
当使用sym函数创建符号变量时，有理算术运算将被启动。
例如，对于双精度矩阵A，sym函数将转化为对应的符号矩阵形式。
A=
1.10001.20001.3000
2.10002.20002.3000
3.10003.20003.3000
A=sym（A）
A=
[11/10,6/5,13/10]
[21/10,11/5,23/10]
[31/10,16/5,33/10]
对于非数值的符号矩阵，通过sym函数，可用与普通矩阵定义类似的格式来定义。
如，下例将定义符号矩阵A、符号矩阵B、符号矩阵C：
Symabcdx
A=sym（‘[ab;cd]’）
A=
[a,b]
[c,d]
B=sym（‘[2*a+b,3*b;5*c+d,2*d]’）
B=
[2*a+b,3*b]
[5*c+d,2*d]
C=sym（‘[1/（1+a）,sin（x）,（b-x）/（a+x）;1,exp（x）,x^2]’）
C=
[1/（1+a）,sin（x）,（b-x）/（a+x）]
[1,exp（x）,x^2]
2．符号矩阵的加、减运算
实现符号矩阵加、减的指令是：
symadd（A,B）%给出两个符号矩阵的和A+B
symsub（A,B）%给出两个符号矩阵的差A-B
显然，参与符号加、减运算的符号矩阵应具有相同的维数。
如，对于上面定义的符号矩阵A和B，利用加、减指令的结果是：
symadd（A,B）
ans=
[3*a+b,4*b]
[6*c+d,3*d]
symsub（A,B）
ans=
[-a-b,-2*b]
[-4*c-d,-d]
对于具有数值的符号矩阵的加法，使用同样计算方法。
如计算两个符号矩阵的和：
A=hilb（3）
A=
1.00000.50000.3333
0.50000.33330.2500
0.33330.25000.2000
A=sym（A）
A=
[1,1/2,1/3]
[1/2,1/3,1/4]
[1/3,1/4,1/5]
B=sym（‘[2*a,b,c;2*b,c;a,b,2*c]’）
B=
[2*a,b,c]
[a,2*b,c]
[a,b,2*c]
symadd（A,B）
ans=
[1+2*a,1/2+b,1/3+c]
[1/2+a,1/3+2*b,1/4+c]
[1/3+a,1/4+b,1/5+2*c]
3．符号矩阵的乘、除运算
实现符号矩阵乘、除的指令是：
symmul（A,B）%给出两个符号矩阵的积A*B
symdiv（A,B）%给出两个符号矩阵的商A/B
如，对于上面定义的符号矩阵A和B，利用乘、除指令的结果是：
symmul（A,B）
ans=
[a*（2*a+b）+b*（5*c+d）,3*a*b+2*b*d]
[c*（2*a+b）+d*（5*c+d）,3*c*b+2*b^2]
symdiv（A,B）
ans=
[（2*a*d-b*d-5*c*d）/（-15*c*b-b*d+4*a*d）,-b*（-b+a）/（-15*c*b-b*d+4*a*d）]
[-d*（3*c+d）/（-15*c*b-b*d+4*a*d）,（2*a*d-3*c*b-b*d）/（-15*c*b-b*d+4*a*d）]
4．符号变量替换
有时，一个符号解或一个符号表达式中，需将一些符号变量替换成数字或其它符号，可利用subs函数实现。
subs函数适用于单个符号矩阵、符号表达式、符号代数方程和微分方程。
subs函数有两种格式：
（1）subs（S,NEW）用新变量NEW替换S中的默认变量。
如：
符号矩阵G，用’pi/3’替换G中的默认变量x：
G=sym（’[a*sin（b+x）,a+b,exp（a*x）,sqrt（x）]’）;
G1=subs（G,’pi/3’）
G1=
[a*sin（b+pi/3）,a+b,exp（a*pi/3）,sqrt（pi/3）]
和利用pi/3替换G中的默认变量x比较：
G=sym（’[a*sin（b+x）,a+b,exp（a*x）,sqrt（x）]’）;
G1=subs（G,pi/3）
G1=
[a*sin（b+1/3*pi）,a+b,exp（1/3*a*pi）,1/3*3^（1/2）*pi^（1/2）]
（2）subs（S,NEW,OLD）用新变量NEW替换S中的指定的变量OLD。
如：
符号表达式f，用2替换f中的变量a：
f=sym（’sin（1/3*a*pi）’）;
subs（f,’2’,’a’）
ans=
sin（1/3*2*pi）
如果用数值替换符号表达式中的变量，将得到相应变量由数值替代的表达式或数值。
如：
f=sym（’sin（1/3*a*pi）’）;
subs（f,2,’a’）
ans=
0.8660
3-4符号微积分运算
1．确定符号变量
当进行数学运算时，对应变量的选取很容易得到。
如，对于函数f=xn，当对f求导时，自然地是对x求导，n看成常数。
而在MATLAB中，如何知道是对x求导而不是对n求导呢？
它通过符号表达式中隐含的符号变量来确定。
在SymbolicMathToolboxs中，确定一个符号表达式中的符号变量的规则是：
（1）只对（除i,j外）单个小写英文字母进行检索。
（2）小写字母x是首选符号变量。
（3）其余小写字母被选用符号变量的次序：
在英文字母表中，靠近“x”的优先，在“x”之后的优先。
按照这一规则，对f=xn求导时，自然是对x求导，n看成常数。
工具箱中还提供了findsym函数来确定表达式中的符号变量。
如，findsym（f,1）寻找第一个符号变量。
findsym函数中的第二个参数，代表了在符号对象中想要寻找的符号变量的个数。
默认时，将给出符号表达式中的所有符号变量。
如下例给出了f中的符号变量。
symsabcx
f=sym（‘a*x^2+b*x+c’）;
findsym（f,1）
ans=
x
findsym（f,2）
ans=
x,c
2．符号微分运算
对符号表达式微分的函数是diff（）。
该函数可以求符号表达式的一阶导数、n阶导数。
调用格式：
（1）diff（f）%传回f对变量x的一次微分值
symsax;
f=sin（a*x）;
df=diff（f）
df=cos（a*x）*a
或：
f=’sin（a*x）’
df=diff（f）
df=
cos（a*x）*a
（2）diff（f,a）%传回f对指定变量a的一次微分值
下列命令分别计算f对变量x和n的微分：
symsxn
f=x^n
diff（f,x）
ans=
x^n*n/x
diff（f,n）
ans=
x^n*log（x）
（3）diff（f,n）%传回f变量x的n次微分值
diff（f,a,n）%传回f对指定变量a的n次微分值
symsax;
f=sin（a*x）;
df=diff（f,x,2）
df=
-sin（a*x）*a^2
diff函数也可以使用符号矩阵作为它的输入，此时，微分按矩阵元素逐个进行。
数值微分函数也是用diff，因此这个函数是靠输入的参数决定是以数值或是符号微分，如果参数为向量则执行数值微分，如果参数为符号表示式则执行符号微分。
先定义下列三个方程式，接着再演算其微分项：
S1='6*x^3-4*x^2+b*x-5';
S2='sin（a）';
S3='（1-t^3）/（1+t^4）';
diff（S1）
ans=
18*x^2-8*x+b
diff（S1,2）
ans=
36*x-8
diff（S1,'b'）
ans=
x
diff（S2）
ans=
cos（a）
diff（S3）
ans=
-3*t^2/（1+t^4）-4*（1-t^3）/（1+t^4）^2*t^3
simplify（diff（S3））
ans=
t^2*（-3+t^4-4*t）/（1+t^4）^2
3．符号积分运算
int函数用以计算函数的积分项，这个函数要找出一符号式F使得diff（F）=f。
如果积
分式的解析式（analyticalform,closedform）不存在的话或是MATLAB无法找到，则int传回原输入的符号式。
相关的函数语法有：
（1）int（f）%传回符号表达式f对变量x的不定积分
symsx;
s=’1/（1+x^2）’;
f=int（f）
f=
atan（x）
（2）int（f,t）%传回符号表达式f对指定变量t的不定积分
symsa;
s=’sin（a*u）’;
int（s,a）
ans=
-1/u*cos（a*u）
（3）int（f,a,b）%传回符号表达式f对变量x的积分值，积分区间为[a,b]，a和b为数值式
int（f,'t',a,b）%传回符号表达式f对变量t的积分值，积分区间为[a,b]，a和b为数值式
int（f,'m','n'）%传回f对变量x的积分值，积分区间为[m,n]，m和n为符号式
我们示范几个例子：
S1='6*x^3-4*x^2+b*x-5';
S2='sin（a）';
S3='sqrt（x）';
int（S1）
ans=
3/2*x^4-4/3*x^3+1/2*b*x^2-5*x
int（S2）
ans=
-cos（a）
int（S3）
ans=
2/3*x^（3/2）
int（S3,'a','b'）
ans=
2/3*b^（3/2）-2/3*a^（3/2）
int（S3,0.5,0.6）
ans=
2/25*15^（1/2）-1/6*2^（1/2）
4．符号微积分运算示例
例1：
对于函数
，求
。
symsx,y;
f=’x*siny’;
dfdxdy=diff（diff（f,x）,y）
dfdxdy=
cos（y）
例2：
对于函数
，先求s关于x的不定积分，再求所得结果关于y的不定积分，即计算
。
symsx,y;
s=’x*exp（-x*y）’;
int（int（s,x）,y）;
ans=
1/y*exp（-x*y）
即
。
例3：
给定一个函数：
，求该函数的微分，然后再将微分结果进行积分运算。
symsxyy1y2
y=’1/（x^2+4*x+3）*sinx’;
y1=diff（y,x）
y1=
-1/（x^2+4*x+3）^2*sinx*（2*x+4）
y2=int（y1,x）
y2=
1/（x^2+4*x+3）*sinx
3-5符号绘图函数ezplot（）
例：
创建符号函数
，并绘该函数的波形。
x=sym（‘x’）;
f=sym（‘1/（5+4*cos（x））’）;
ezplot（f）
求f的微分并绘出结果的波形
f1=diff（f）
f1=
4/（5+4*cos（x））^2*sin（x）
ezplot（f1）
求f的二次微分并绘出结果的波形
f2=diff（f,2）
f2=
32/（5+4*cos（x））^3*sin（x）^2+4/（5+4*cos（x））^2*cos（x）
ezplot（f2）
求f2的二次积分并绘出结果的波形
g=int（int（f2））
g=
-8/（tan（1/2*x）^2+9）
ezplot（g）
此时，我们发现，g和原函数f形式不同，但图形相似。
下面求g和f的差值：
r=f-g
r=
1/（5+4*cos（x））+8/（tan（1/2*x）^2+9）
ezplot（r）
可见，图形窗口中绘制的f-g的波形是一条幅度为1的直线。
对r进行简化，
e=simple（r）
e=
1
所以，f-g为是个恒定常数。
3-6求解常微分方程式
MATLAB解常微分方程式的语法是dsolve（'equation','condition'），其中equation代表常微分方程式即y'=g（x,y），且须以Dy代表一阶微分项y' D2y代表二阶微分项y''，condition则为初始条件。
假设有以下三个一阶常微分方程式和其初始条件
y'=3x2，y
（2）=0.5
y'=2xcos（y）2，y（0）=0.25
y'=3y+exp（2x），y（0）=3
对应上述常微分方程式的符号运算式为：
soln_1=dsolve（'Dy=3*x^2','y
（2）=0.5'）
ans=
x^3-7.500000000000000
ezplot（soln_1,[2,4]）%看看这个函数的长相
soln_2=dsolve（'Dy=2*x*cos（y）^2','y（0）=pi/4'）
ans=
atan（x^2+1）
soln_3=dsolve（'Dy=3*y+exp（2*x）','y（0）=3'）
ans=
-exp（2*x）+4*exp（3*x）
3-7连续信号的微分与积分
连续信号的微分可由diff近似计算：
h=0.001;x=0:
h:
pi;
y=diff（sin（x.^2））/h;
定积分可由quad函数或quadl函数实现:
quad（‘function_name’,a,b）%function_name为被积函数名，a和b指积分区间
例：
求example2-5（第二讲中）所示三角波求微分和积分
[第二讲中三角波的产生：
%example2-5
t=-4:
0.001:
4;
yt=tripuls（t,4,0.5）;
plot（t,yt）]
先将三角波写成MATLAB函数（m文件），名为x2_5:
functionyt=x2_5（t）
yt=tripuls（t,4,0.5）;
用diff和quad函数可实现三角波的微分和积分
%example3-1differentiation（微分）
h=0.001;x=-3:
h:
3;
y1=diff（x2_5（t））*1/h;
plot（t（1:
length（t）-1）,y1）
title（‘dx（t）/dt’））
%example3-2integration（积分）
t=-3:
0.1:
3;
forn=1:
length（t）
y2（n）=quad（'x1_2',-3,t（n））;
end
plot（t,y2）
title（'integralofx（t）'）
3-8练习
1、用符号表达式表示信号
，并绘制波形。
2、实现下列连续时间信号x（t）的微分和积分。
3、信号
，用MATLAB符号运算的相关命令求：
、
、
、
、
，并绘出其时域波形。
[3答案：
symst%定义符号变量t
f=sym（’（t/2+1）*（heaviside（t+2）-heaviside（t-2））’）%创建符号表达式f，heaviside函数见第二讲
subplot（2,3,1）,ezplot（f,[-3,3]）%subplot把图形窗口分成2行3列的6个小窗口，当前波形画在第1个小窗口中，下同
y1=subs（f,t,t+2）%subs（）函数指将连续信号f中的时间变量t用t+2替换，下同
subplot（2,3,2）,ezplot（y1,[-5,1]）
y2=subs（f,t,t-2）
subplot（2,3,3）,ezplot（y2,[-1,5]）
y3=subs（f,t,-t）
subplot（2,3,4）,ezplot（y3,[-3,3]）
y4=subs（f,t,2*t）
subplot（2,3,5）,ezplot（y4,[-2,2]）
y5=-f
subplot（2,3,6）,ezplot（y5,[-3,3]）
请注意：
上述命令执行后，f、y1、y2、y3、y4、y5的符号表达式（在Command窗口中显示）。
]
升级会员 