第3讲 matlab的符号运算Word文档格式.docx
《第3讲 matlab的符号运算Word文档格式.docx》由会员分享,可在线阅读,更多相关《第3讲 matlab的符号运算Word文档格式.docx(21页珍藏版)》请在冰点文库上搜索。
f=sym(‘a*x^2+b*x+c’)
f=‘a*x^2+b*x+c’
此例中,将符号表达式赋给符号变量f,但这不是必需的,引入符号变量是为了以后调用方便。
在这种情况下,没有创建对应于表达式中a、b、c、x项的变量,为了执行符号数学运算(如微分、积分等),必须显式地创建这些变量,可用下列命令创建:
symsabcx
如下例中创建了符号表达式和符号方程,分别赋给相应的符号对象。
symsxabc
%创建符号表达式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’)
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)
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)
[1,1/2,1/3]
[1/2,1/3,1/4]
[1/3,1/4,1/5]
当调用SymbolicMathToolbox时,默认地调用有理算术运算。
当使用sym函数创建符号变量时,有理算术运算将被启动。
例如,对于双精度矩阵A,sym函数将转化为对应的符号矩阵形式。
1.10001.20001.3000
2.10002.20002.3000
3.10003.20003.3000
[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,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)
[3*a+b,4*b]
[6*c+d,3*d]
symsub(A,B)
[-a-b,-2*b]
[-4*c-d,-d]
对于具有数值的符号矩阵的加法,使用同样计算方法。
如计算两个符号矩阵的和:
[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]’)
[2*a,b,c]
[a,2*b,c]
[a,b,2*c]
[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)
[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)
[(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比较:
G1=subs(G,pi/3)
[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’)
sin(1/3*2*pi)
如果用数值替换符号表达式中的变量,将得到相应变量由数值替代的表达式或数值。
subs(f,2,’a’)
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中的符号变量。
f=sym(‘a*x^2+b*x+c’);
findsym(f,1)
x
findsym(f,2)
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=
cos(a*x)*a
(2)diff(f,a)%传回f对指定变量a的一次微分值
下列命令分别计算f对变量x和n的微分:
symsxn
f=x^n
diff(f,x)
x^n*n/x
diff(f,n)
x^n*log(x)
(3)diff(f,n)%传回f变量x的n次微分值
diff(f,a,n)%传回f对指定变量a的n次微分值
df=diff(f,x,2)
-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)
18*x^2-8*x+b
diff(S1,2)
ans=
36*x-8
diff(S1,'
b'
)
x
diff(S2)
cos(a)
diff(S3)
-3*t^2/(1+t^4)-4*(1-t^3)/(1+t^4)^2*t^3
simplify(diff(S3))
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)
-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为数值式
m'
'
n'
)%传回f对变量x的积分值,积分区间为[m,n],m和n为符号式
我们示范几个例子:
sqrt(x)'
int(S1)
3/2*x^4-4/3*x^3+1/2*b*x^2-5*x
int(S2)
-cos(a)
int(S3)
2/3*x^(3/2)
int(S3,'
a'
2/3*b^(3/2)-2/3*a^(3/2)
int(S3,0.5,0.6)
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的不定积分,即计算
s=’x*exp(-x*y)’;
int(int(s,x),y);
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
=2xcos(y)2,y(0)=0.25
=3y+exp(2x),y(0)=3
对应上述常微分方程式的符号运算式为:
soln_1=dsolve('
Dy=3*x^2'
y
(2)=0.5'
x^3-7.500000000000000
ezplot(soln_1,[2,4])%看看这个函数的长相
soln_2=dsolve('
Dy=2*x*cos(y)^2'
y(0)=pi/4'
atan(x^2+1)
soln_3=dsolve('
Dy=3*y+exp(2*x)'
y(0)=3'
-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)
用diff和quad函数可实现三角波的微分和积分
%example3-1differentiation(微分)
x=-3:
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:
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窗口中显示)。
]