Matlab 参考教程第六章 符号计算.docx
- 文档编号:8998559
- 上传时间:2023-02-02
- 格式:DOCX
- 页数:30
- 大小:269.79KB
Matlab 参考教程第六章 符号计算.docx
《Matlab 参考教程第六章 符号计算.docx》由会员分享,可在线阅读,更多相关《Matlab 参考教程第六章 符号计算.docx(30页珍藏版)》请在冰豆网上搜索。
Matlab参考教程第六章符号计算
6符号计算
MATLAB4.2中,符号计算所依赖的SymbolicMathToolbox1.0版是一个过渡性产品。
1.0版中的几乎所有指令都已经被废止。
而今MATLAB5.3的符号计算工具包已升级为2.1版,它的工作原动机是MapleV5。
2.1版采用全新的数据结构、面向对象编程和重载技术,使得符号计算和数值计算在形式和风格上浑然统一。
6.1符号对象和符号表达式
6.1.1符号对象的生成和使用
【*例6.1.1-1】符号常数形成中的差异
a1=[1/3,pi/7,sqrt(5),pi+sqrt(5)]%a1是数值常数<1>
a2=sym([1/3,pi/7,sqrt(5),pi+sqrt(5)])%最接近的有理表示<2>
a3=sym([1/3,pi/7,sqrt(5),pi+sqrt(5)],'e')%带估计误差的有理表示<3>
a4=sym('[1/3,pi/7,sqrt(5),pi+sqrt(5)]')%绝对准确的符号数值表示<4>
a24=a2-a4
a1=
0.33330.44882.23615.3777
a2=
[1/3,pi/7,sqrt(5),6054707603575008*2^(-50)]
a3=
[1/3-eps/12,pi/7-13*eps/165,sqrt(5)+137*eps/280,6054707603575008*2^(-50)]
a4=
[1/3,pi/7,sqrt(5),pi+sqrt(5)]
a24=
[0,0,0,189209612611719/35184372088832-pi-5^(1/2)]
【*例6.1.1-2】演示:
几种输入下产生矩阵的异同。
a1=sym([1/3,0.2+sqrt
(2),pi])%产生
符号数组<1>
a2=sym('[1/3,0.2+sqrt
(2),pi]')%产生
符号数组<2>
a3=sym('[1/30.2+sqrt
(2)pi]')%2.1版中产生
符号数组<3>
a1_a2=a1-a2%为比较a1,a2
a1=
[1/3,7269771597999872*2^(-52),pi]
a2=
[1/3,0.2+sqrt
(2),pi]
a3=
[1/3,0.2+sqrt
(2)pi]
a1_a2=
[0,1.4142135623730951010657008737326-2^(1/2),0]
【*例6.1.1-3】把字符表达式转换为符号变量
y=sym('2*sin(x)*cos(x)')%把字符表达式转换为符号变量
y=simple(y)%按规则把已有的y符号表达式化成最简形式
y=
2*sin(x)*cos(x)
y=
sin(2*x)
【*例6.1.1-4】用符号计算验证三角等式
。
symsfai1fai2;y=simple(sin(fai1)*cos(fai2)-cos(fai1)*sin(fai2))
y=
sin(fai1-fai2)
【*例6.1.1-5】求矩阵
的行列式值、逆和特征根
symsa11a12a21a22;A=[a11,a12;a21,a22]
DA=det(A),IA=inv(A),EA=eig(A)
A=
[a11,a12]
[a21,a22]
DA=
a11*a22-a12*a21
IA=
[a22/(a11*a22-a12*a21),-a12/(a11*a22-a12*a21)]
[-a21/(a11*a22-a12*a21),a11/(a11*a22-a12*a21)]
EA=
[1/2*a11+1/2*a22+1/2*(a11^2-2*a11*a22+a22^2+4*a12*a21)^(1/2)]
[1/2*a11+1/2*a22-1/2*(a11^2-2*a11*a22+a22^2+4*a12*a21)^(1/2)]
【*例6.1.1-6】验证积分
。
symsAttaow;yf=int(A*exp(-i*w*t),t,-tao/2,tao/2);Yf=simple(yf)
Yf=
2*A*sin(1/2*tao*w)/w
6.1.2符号计算中的算符和基本函数
6.1.3识别对象类别的指令
【*例6.1.3-1】数据对象及其识别指令的使用。
(1)生成三种不同类型的矩阵,给出不同的显示形式
clear,a=1;b=2;c=3;d=4;%产生四个数值变量
Mn=[a,b;c,d]%利用已赋值变量构成数值矩阵
Mc='[a,b;c,d]'%字符串中的a,b,c,d与前面输入的数值变量无关
Ms=sym(Mc)%Ms是一个符号变量,它与前面各变量无关。
Mn=
12
34
Mc=
[a,b;c,d]
Ms=
[a,b]
[c,d]
(2)三种矩阵的大小不同
SizeMn=size(Mn),SizeMc=size(Mc),SizeMs=size(Ms)
SizeMn=
22
SizeMc=
19
SizeMs=
22
(3)用class获得每种矩阵的类别
CMn=class(Mn),CMc=class(Mc),CMs=class(Ms)
CMn=
double
CMc=
char
CMs=
sym
(4)用isa判断每种矩阵的类别(若返回1,表示判断正确)
isa(Mn,'double'),isa(Mc,'char'),isa(Ms,'sym')
ans=
1
ans=
1
ans=
1
(5)利用whos观察内存变量的类别和其它属性
whosMnMcMs%观察三个变量的类别和属性
NameSizeBytesClass
Mc1x918chararray
Mn2x232doublearray
Ms2x2408symobject
Grandtotalis21elementsusing458bytes
6.1.4符号表达式中自由变量的确定
【*例6.1.4-1】对独立自由符号变量的自动辨认。
(1)生成符号变量
symsabxXY;k=sym('3');z=sym('c*sqrt(delta)+y*sin(theta)');
EXPR=a*z*X+(b*x^2+k)*Y;
(2)找出EXPR中的全部自由符号变量
findsym(EXPR)%除常数符号k外的所有独立符号变量都被列出
ans=
X,Y,a,b,c,delta,theta,x,y
(3)在EXPR中确定一个自由符号变量
findsym(EXPR,1)
ans=
x
(4)在EXPR中确定2个和3个自由变量时的执行情况
findsym(EXPR,2),findsym(EXPR,3)
ans=
x,y
ans=
x,y,theta
【*例6.1.4-2】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]
ans=
x
6.2符号表达式和符号函数的操作
6.2.1符号表达式的操作
【*例6.2.1-1】按不同的方式合并同幂项。
EXPR=sym('(x^2+x*exp(-t)+1)*(x+exp(-t))');
expr1=collect(EXPR)%默认合并x同幂项系数
expr2=collect(EXPR,'exp(-t)')%合并exp(-t)同幂项系数
expr1=
x^3+2*exp(-t)*x^2+(1+exp(-t)^2)*x+exp(-t)
expr2=
x*exp(-t)^2+(2*x^2+1)*exp(-t)+(x^2+1)*x
【*例6.2.1-2】factor指令的使用
(1)除x外不含其它自由变量的情况
symsax;f1=x^4-5*x^3+5*x^2+5*x-6;factor(f1)
ans=
(x-1)*(x-2)*(x-3)*(x+1)
(2)含其它自由变量的情况之一
f2=x^2-a^2;factor(f2)
ans=
(x-a)*(x+a)
(3)对正整数的质数分解
factor(1025)
ans=
5541
【*例6.2.1-3】对多项式进行嵌套型分解
clear;symsax;f1=x^4-5*x^3+5*x^2+5*x-6;horner(f1)
ans=
-6+(5+(5+(-5+x)*x)*x)*x
【*例6.2.1-4】写出矩阵
各元素的分子、分母多项式
(1)求矩阵各元素的分子、分母多项式
symsx;A=[3/2,(x^2+3)/(2*x-1)+3*x/(x-1);4/x^2,3*x+4];
[n,d]=numden(A)
pretty(simplify(A))%为方便阅读而设。
请用simple代替simplify试试。
<3>
n=
[3,x^3+5*x^2-3]
[4,3*x+4]
d=
[2,(2*x-1)*(x-1)]
[x^2,1]
[32]
[x+5x-3]
[3/2-----------------]
[(2x-1)(x-1)]
[]
[4]
[----3x+4]
[2]
[x]
(2)请读者自己运行以下指令,显示结果应与指令<3>相同。
pretty(simplify(n./d))%注意这里用的是“数组除”,而不是“矩阵除”
[32]
[x+5x-3]
[3/2-----------------]
[(2x-1)(x-1)]
[]
[4]
[----3x+4]
[2]
[x]
【*例6.2.1-5】简化
(1)运用simplify简化(即使多次运用simplify也不能得到最简形式。
)
symsx;f=(1/x^3+6/x^2+12/x+8)^(1/3);
sfy1=simplify(f),sfy2=simplify(sfy1)
sfy1=
((2*x+1)^3/x^3)^(1/3)
sfy2=
((2*x+1)^3/x^3)^(1/3)
(2)运用simple简化
g1=simple(f),g2=simple(g1)
g1=
(2*x+1)/x
g2=
2+1/x
【*例6.2.1-6】简化
symsx;ff=cos(x)+sqrt(-sin(x)^2);
ssfy1=simplify(ff),ssfy2=simplify(ssfy1)
ssfy1=
cos(x)+(-sin(x)^2)^(1/2)
ssfy2=
cos(x)+(-sin(x)^2)^(1/2)
gg1=simple(ff),gg2=simple(gg1)
gg1=
cos(x)+i*sin(x)
gg2=
exp(i*x)
6.2.2符号函数的求反和复合
【*例6.2.2-1】求
的反函数
symsx;f=x^2;g=finverse(f)
Warning:
finverse(x^2)isnotunique.
>InE:
\MAT53\toolbox\symbolic\@sym\finverse.matline43
g=
x^(1/2)
fg=simple(compose(g,f))%验算g(f(x))是否等于x
fg=
x
【*例6.2.2-2】求
的复合函数
(1)自变量由机器确定
symsxyufait;f=x/(1+u^2);g=cos(y+fai);fg1=compose(f,g)
fg1=
cos(y+fai)/(1+u^2)
(2)指定自变量
fg2=compose(f,g,u,fai,t)
fg2=
x/(1+cos(y+t)^2)
6.2.3置换及其应用
6.2.3.1自动执行的子表达式置换指令
【*例6.2.3.1-1】演示子表达式的置换表示。
clearall,symsabcdW;[V,D]=eig([ab;cd]);
[RVD,W]=subexpr([V;D],W)%<2>
RVD=
[-(1/2*d-1/2*a-1/2*W)/c,-(1/2*d-1/2*a+1/2*W)/c]
[1,1]
[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)
6.2.3.2通用置换指令
【*例6.2.3.2-1】用简单算例演示subs的置换规则。
(1)产生符号函数
symsax;f=a*sin(x)+5;
(2)符号变量置换
f1=subs(f,'sin(x)',sym('y'))%<2>
f1=
a*y+5
(3)符号常数置换
f2=subs(f,{a,x},{2,sym(pi/3)})%<3>
f2=
3^(1/2)+5
(4)双精度数值置换(即所有自由变量被双精度数值取代。
取a=2,x=pi/3)
f3=subs(f,{a,x},{2,pi/3})%<4>
f3=
6.7321
(5)数值数组置换之一(取a=2,x=0:
pi/6:
pi)
f4=subs(subs(f,a,2),x,0:
pi/6:
pi)%<5>
f4=
5.00006.00006.73217.00006.73216.00005.0000
(6)数值数组置换之二(取a=0:
6,x=0:
pi/6:
pi)
f5=subs(f,{a,x},{0:
6,0:
pi/6:
pi})%<6>
f5=
5.00005.50006.73218.00008.46417.50005.0000
6.2.4符号数值精度控制和任意精度计算
6.2.4.1向双精度数值转换的doblue指令
6.2.4.2任意精度的符号数值
【*例6.2.4.2-1】指令使用演示。
digits%显示省缺符号数值计算相对精度
Digits=32
p0=sym('(1+sqrt(5))/2');%p0为(1+sqrt(5))/2准确值
p1=sym((1+sqrt(5))/2)%p1是(1+sqrt(5))/2在数值环境下的近似值
e01=vpa(abs(p0-p1))%在符号环境32位省缺精度下,观察p0,p1间误差
p1=
7286977268806824*2^(-52)
e01=
.543211520368250e-16
p2=vpa(p0)%p2是p0在32位相对精度下的近似
e02=vpa(abs(p0-p2),64)%在64位相对精度下,观察p0,p2间的误差
p2=
1.6180339887498948482045868343657
e02=
.61882279690820194237137864551377e-31
digits%验证vpa变精度计算不影响全局的符号数值计算精度
Digits=32
6.2.5符号对象与其它数据对象间的转换
【*例6.2.5-1】符号、数值间的转换。
phi=sym((1+sqrt(5))/2)%把数值对象转换成符号常数
double(phi)%把符号常数转换为双精度存储的数值
phi=
7286977268806824*2^(-52)
ans=
1.6180
【*例6.2.5-2】各种多项式表示形式之间的转换
symsx;f=x^3+2*x^2-3*x+5;%生成符号多项式
sy2p=sym2poly(f)%由符号多项式产生数值系数行向量
p2st=poly2str(sy2p,'x')%把系数行向量变成易读表示式
p2sy=poly2sym(sy2p)%把数值系数行向量再转换为符号多项式
pretty(f,'x')%显示符号多项式的易读表示形式
sy2p=
12-35
p2st=
x^3+2x^2-3x+5
p2sy=
x^3+2*x^2-3*x+5
32
x+2x-3x+5
6.3符号微积分
6.3.1符号序列的求和
【*例6.3.1-1】求
,
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))%f2的自变量被确认为k
s1=
[1/2*t*(t-1),k^3*t]
s2=
[1/8*pi^2,-log
(2)]
6.3.2符号微分和
矩阵
【*例6.3.2-1】求
、
和
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]
[0,0]
dfdxdt=
[0,0]
[-sin(x),0]
【*例6.3.2-2】求
的
矩阵。
symsx1x2x3;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)]
6.3.3符号积分
6.3.3.1通用积分指令
6.3.3.2交互式近似积分指令
6.3.3.3符号积分示例
【*例6.3.3.3-1】求
。
演示:
积分指令对符号函数矩阵的作用。
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)]
【*例6.3.3.3-2】求
。
演示如何使用mfun指令获取一组积分值。
(1)求一般积分结果
F1=int('1/log(t)','t',0,'x')
F1=
-Ei(1,-log(x))
(2)利用mfun指令求x=0.5,0.6,0.7,0.8,0.9时的定积分
x=0.5:
0.1:
0.9
F115=-mfun('Ei',1,-log(x))
x=
0.50000.60000.70000.80000.9000
F115=
-0.3787-0.5469-0.7809-1.1340-1.7758
【*例6.3.3.3-3】求积分
。
注意:
内积分上下限都是函数。
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-6072064/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2^(3/4)
VF2=
224.92153573331143159790710032805
【*例6.3.3.3-4】利用rsums求
积分。
(与例6.3.3.3-2结果比较)
symsxpositive;px=0.5/log(0.5*x);rsums(px)
图6.3.3.3-4交互式近似积分
6.3.4符号卷积
【*例6.3.4-1】本例演示卷积的时域积分法:
已知系统冲激响应
,求
输入下的输出响应。
symsTttao;ut=exp(-t);%定义系统输入
ht=exp(-t/T)/T;%定义系统冲激响应
uh_tao=subs(ut,t,tao)*subs(ht,t,t-tao);%运用变量替换指令形成被积函数
yt=int(uh_tao,tao,0,t);%实施卷积
yt=simple(yt)
yt=
-(exp(-t)-exp(-t/T))/(T-1)
【*例6.3.4-2】本例演示通过变换和反变换求取卷积。
系统冲激响应、输入同上例,求输出。
symss;yt=ilaplace(laplace(ut,t,s)*laplace(ht,t,s),s,t);yt=simple(yt)
yt=
-(exp(-t)-exp(-t/T))/(T-1)
【*例6.3.4-3】求函数
和
的卷积。
(1)在5.2版(配SymbolicMathToolbox2.0.1)中,采用以下指令。
symsttao;ut=sym('Heaviside(t)-Heaviside(t-1)');ht=t*exp(-t);
yt=int(subs(ut,t,tao)*subs(ht,t,abs(t)-tao),tao,0,abs(t));
yt=collect(yt,'signum(abs(t)-1)'),yt=subs(yt,abs(t),t)%<3>
yt
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 参考教程第六章 符号计算 参考 教程 第六 符号 计算