ch05符号计算总结.docx
《ch05符号计算总结.docx》由会员分享,可在线阅读,更多相关《ch05符号计算总结.docx(93页珍藏版)》请在冰豆网上搜索。
ch05符号计算总结
第5章符号计算
符号计算的优点:
凭借恒等式,数学定理,通过推理和演绎,给出具有“无限尺度”描写能力的解析结果。
当没有封闭解时,符号计算则妥协地给出任意精度数值解。
所谓“无限尺度”是相对数值计算的“有限精度、有限空间”而言的。
MATLAB本身由数值计算引擎驱动,而没有符号演绎能力。
2008年前,MATLAB的符号计算能力借助于Maple。
现在,随MATLAB默认安装的符号计算引擎是MuPAD。
虽然MATLAB尽力地保持着形式上的向前兼容,但引擎换装确实导致MATLAB符号计算环境发生了根本改变。
本章内容完全针对MuPAD引擎的符号计算而写,分三个层面:
●无需任何MuPAD知识,仅使用MATLAB符号数学工具包(Symbolicmathtoolbox)提供的(前台)函数实施的符号计算及仿真。
本章第1节比较完整地描述了符号计算的机理、规则和帮助系统。
第2到第6节所涉内容包括微积分、微分方程、积分变换、矩阵分析和代数方程。
第7节代数状态方程和第8节数据探索的内容,用以表现现代计算能力对传统方法或技巧的冲击。
●需要少量MuPAD知识,拓展符号计算应用范围的内容安排在第5.9节。
在这一节中不仅讲述特殊函数计算,而且还介绍了如何利用符号函数产生M函数文件,如何利用符号函数制作用户所需的Simulink模块。
●需要较多MuPAD知识,借助evalin和feval指令进入MuPAD空间完成符号计算。
对于那些在MATLAB中仍选择Maple符号计算引擎的读者来说,本章第一层面的内容仍可借鉴。
值得指出:
随书光盘mbook目录上的“ch05_符号计算.doc”保存有本章全部算例的运作指令、计算结果和图形;mfile目录上则保存有本章算例中所有带exm前缀文件名的M文件和MEX文件。
.1符号对象的产生和识别
.1.1基本符号对象的创建
101定义符号数字和符号常数
102定义基本符号变量
103定义元符号表达式
.1.2符号计算中的算符和函数指令
101符号计算中的算符
102符号计算中的函数指令
.1.3符号对象、变量、自由变量的识别
101符号对象的识别
【例5.1-1】数据对象及其识别指令的使用。
(1)
clear
a=1;b=2;c=3;d=4;
Mn=[a,b;c,d]
Mc='[a,b;c,d]'
Ms=sym(Mc)
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)
CMn=class(Mn)
CMc=class(Mc)
CMs=class(Ms)
CMn=
double
CMc=
char
CMs=
sym
(4)
isa(Mn,'double')
isa(Mc,'char')
isa(Ms,'sym')
ans=
1
ans=
1
ans=
1
(5)
whosMnMcMs
NameSizeBytesClassAttributes
Mc1x918char
Mn2x232double
Ms2x260sym
102符号变量及自由变量的认定
【例5.1-2】基本符号变量和衍生符号表达式的定义、符号变量、自由符号变量的机器辨认。
(1)
clear
symsabktwxyzXA
c=sym('22');
f=sym('M*N');
ex1=c+f+(exp(-a*t)*sin(w*t)+b*y+k*x+A*X+log(z))
ex1=
log(z)+b*y+k*x+sin(t*w)/exp(a*t)+A*X+M*N+22
(2)
symvar(ex1)%<6>
ans=
[A,M,N,X,a,b,k,t,w,x,y,z]
findsym(ex1)
ans=
A,M,N,X,a,b,k,t,w,x,y,z
(3)
symvar(ex1,20)%<8>
ans=
[x,y,w,z,t,k,b,a,X,N,M,A]
findsym(ex1,13)%
ans=
x,y,w,z,t,k,b,a,X,N,M,A
(4)
who%<10>
%
Yourvariablesare:
AXaansbcex1fktwxyz
(5)
dex1dx=diff(ex1,z)%<11>
iex1db=int(ex1,b)%<12>
dex1dx=
1/z
iex1db=
(y*b^2)/2+(log(z)+k*x+sin(t*w)/exp(a*t)+A*X+M*N+22)*b
【例5.1-3】元符号表达式的定义、符号变量、自由符号变量的机器辨认。
(1)
clear
c=sym('22');
f=sym('M*N');
ex2=c+f+sym('exp(-a*t)*sin(w*t)+b*y+k*x+A*X+log(z)')
ex2=
log(z)+b*y+k*x+sin(t*w)/exp(a*t)+A*X+M*N+22
(2)
symvar(ex2)%<5>
ans=
[A,M,N,X,a,b,k,t,w,x,y,z]
(3)
symvar(ex2,20)%<6>
ans=
[x,y,w,z,t,k,b,a,X,N,M,A]
(4)
who%<7>
Yourvariablesare:
anscex2f
(5)
dex1dx=diff(ex2,z)%<8>
?
?
?
Undefinedfunctionorvariable'z'.
dex1db=int(ex2,b)%<9>
?
?
?
Undefinedfunctionorvariable'b'.
(6)
E3=sym('a*sqrt(theta)')%
?
?
?
Errorusing==>sym.sym>convertExpressionat2547
Error:
argumentmustbeof'Type:
:
Arithmetical'[sqrt]
Errorin==>sym.sym>convertCharat2458
s=convertExpression(x);
Errorin==>sym.sym>convertCharWithOptionat2441
s=convertChar(x);
Errorin==>sym.sym>tomupadat2195
S=convertCharWithOption(x,a);
Errorin==>sym.sym>sym.symat111
S.s=tomupad(x,'');
E4=sym('a*sqrt(theta123)')
E4=
a*theta123^(1/2)
【例5.1-4】symvar确定自由变量是对整个数组进行的。
(1)
clear
symsabtvwxyz
A=[a+b*x,y*sin(t)+w;x*exp(-t),log(z)+v]
A=
[a+b*x,w+y*sin(t)]
[x/exp(t),v+log(z)]
(2)
symvar(A,1)
ans=
x
symvar(A,3)
ans=
[x,y,w]
symvar(A,10)
ans=
[x,y,w,z,v,t,b,a]
.1.4符号运算机理和变量假设
101符号运算的工作机理
102对符号变量的限定性假设
(1)
(2)
103清除变量和撤销假设
【例5.1-5】syms对变量所作限定性假设的影响。
(1)
symsxclear
f=x^3+4.75*x+2.5;
rf=solve(f,x)
rf=
-1/2
1/4-(79^(1/2)*i)/4
1/4+(79^(1/2)*i)/4
evalin(symengine,'getprop(x)')%<4>
ans=
C_
evalin(symengine,'property:
:
showprops(x);')
%<5>
ans=
[emptysym]
(2)
symsxreal
rfr=solve(f,x)
rfr=
-1/2
evalin(symengine,'getprop(x)')%<8>
ans=
R_
evalin(symengine,'property:
:
showprops(x);')%<9>
ans=
[xinR_]
(3)
clearx
symsx
g=x^2+x+5;
rg=solve(g,x)
Warning:
Explicitsolutioncouldnotbefound.
>Insolveat81
rg=
[emptysym]
(4)
symsxclear
rg=solve(g,x)
rg=
-1/2-(19^(1/2)*i)/2
-1/2+(19^(1/2)*i)/2
【例5.1-6】借助evalin对符号变量进行假设的设定和撤销。
(1)
clear
symsx
f=x^3-5.25*x-2.5;
rf=feval(symengine,'solve',f,'x')%<4>
rf=
[-1/2,5/2,-2]
(2)
x=sym(x,'positive');
rf=feval(symengine,'solve',f,'x')%<6>
rf=
5/2
rfs=solve(f,x)
rfs=
5/2
evalin(symengine,'getprop(x)')%<8>
ans=
Dom:
:
Interval(0,Inf)
evalin(symengine,'deletex')
evalin(symengine,'getprop(x)')%<10>
ans=
C_
(3)
evalin(symengine,'assume(x,Type:
:
Integer);')%<11>
rf=feval(symengine,'solve',f,'x')%<12>
rf=
-2
evalin(symengine,'assumeAlso(x>0)')%<13>
rf_ev=evalin(symengine,'solve(x^3-5.25*x-2.5);')%<14>
rf_evc=evalin(symengine,['solve(',char(f),')'])%<15>
rf_f=feval(symengine,'solve',f,'x')%<16>
rf_ev=
[emptysym]
rf_evc=
[emptysym]
rf_f=
[emptysym]
evalin(symengine,'assume(x,Type:
:
Complex);')%<17>
rf=feval(symengine,'solve',f,'x')%<18>
rf=
[-1/2,5/2,-2]
(4)
syms%<19>
'ans''f''rf''rf_ev''rf_evc''rf_f''rfs''x'
.1.5符号帮助及其他常用指令
101符号运作的帮助体系
(1)“直接调用符号计算指令”的求助
(2)借助mfun调用的“特殊函数指令”的求助
(3)借助evalin或feval调用的“MuPAD指令”的求助。
【例5.1-7】关于laplace,erfc,rec三个指令的求助过程。
(1)
图5.1-1“可直接调用的符号计算指令”的帮助信息搜索
(2)
图5.1-2MATLAB帮助浏览器展示全部特殊函数符号计算指令及用法
(3)
图5.1-3MuPAD帮助浏览器展示的rec指令帮助信息
102服务于符号运作的其他指令
【例5.1-8】各种帮助指令、符号变量罗列指令
(1)
clearall
reset(symengine)%<2>
Da=1.2;Dw=1/3;%
symssaswsxsysz%
symsABpositive%<5>
symsCreal%<6>
(2)
whos%<7>
NameSizeBytesClassAttributes
A1x160sym
B1x160sym
C1x160sym
Da1x18double
Dw1x18double
sa1x160sym
sw1x160sym
sx1x160sym
sy1x160sym
sz1x160sym
(3)
syms%<8>
'A''B''C''sa''sw''sx''sy''sz'
(4)
evalin(symengine,'anames(Properties)')%<9>
ans=
[A,B,C]
(5)
clearA
syms%<11>
'B''C''ans''sa''sw''sx''sy''sz'
evalin(symengine,'anames(Properties)')%<12>
ans=
[A,B,C]
(6)
symsBclear%<13>
syms
'B''C''ans''sa''sw''sx''sy''sz'
evalin(symengine,'anames(Properties)')%<15>
ans=
[A,C]
.2数字类型转换及符号表达式操作
.2.1数字类型及转换
101三种数字类型及转换指令
图5.2-1符号、字符、数值间的相互转换
102双精度数字向符号数字转换
【例5.2-1】本例目的:
准确符号数字的产生;双精度数字转换成符号数字的各种格式;由双精度数字转换而得符号数的误差;vpa的用法;数字类别的判断。
(1)
clear
sa=sym('1/3')%<2>
sb=sym('pi+sqrt(5)')%<3>
sa=
1/3
sb=
pi+5^(1/2)
(2)
formatlong
a=1/3%<5>
b=pi+sqrt(5)%<6>
sa_a=vpa(sa-a)%<7>
sb_b=vpa(sb-b)%<8>
a=
0.333333333333333
b=
5.377660631089583
sa_a=
0.0
sb_b=
0.000000000000000013822375841085200048593542564188
(3)
asr=sym(a)
bsr=sym(b)
sa_asr=vpa(sa-asr)
sb_bsr=vpa(sb-bsr)
asr=
1/3
bsr=
189********1719/35184372088832
sa_asr=
0.0
sb_bsr=
0.000000000000000013822375841085200048593542564188
(4)
ase=sym(a,'e')
bse=sym(b,'e')
sa_ase=vpa(sa-ase)
sb_bse=vpa(sb-bse)
ase=
1/3-eps/12
bse=
189********1719/35184372088832
sa_ase=
0.083333333333333333333333333333333*eps
sb_bse=
0.000000000000000013822375841085200048593542564188
(5)
asf=sym(a,'f')
bsf=sym(b,'f')
sa_asf=vpa(sa-asf)
sb_bsf=vpa(sb-bsf)
asf=
6004799503160661/180********481984
bsf=
189********1719/35184372088832
sa_asf=
0.000000000000000018503717077085942340393861134847
sb_bsf=
0.000000000000000013822375841085200048593542564188
(6)
asd=sym(a,'d')
bsd=sym(b,'d')
sa_asd=vpa(sa-asd)
sb_bsd=vpa(sb-bsd)
asd=
0.33333333333333331482961625624739
bsd=
5.3776606310895829210494412109256
sa_asd=
0.000000000000000018503717077085943333333327037792
sb_bsd=
0.000000000000000013822375841085179119638453173161
(7)
class(sa)
ans=
sym
disp('a的类别asr的类别sa_a的类别')
disp([class(a),blanks(4),class(asr),blanks(7),class(sa_a)])
a的类别asr的类别sa_a的类别
doublesymsym
【例5.2-2】本例目的:
揭露字符串数字在不同计算引擎中的不同表现;指令format的用法;double与str2double的区别。
(1)
clear
formatlong%<2>
ad=1/2+3^(2/3)
astr='1/2+3^(2/3)'
asym=sym(astr)
ad=
2.580083823051904
astr=
1/2+3^(2/3)
asym=
3^(2/3)+1/2
(2)
disp('adastrasym')
disp([blanks(6),'的数据类别'])
disp([class(ad),blanks(3),class(astr),blanks(4),class(asym)])
adastrasym
的数据类别
doublecharsym
(3)
asymPLUSad=asym+sym(0.1)%<9>
asymPLUSda2=asym+0.1%<10>
asymPLUSad=
3^(2/3)+3/5
asymPLUSda2=
3^(2/3)+3/5
(4)
astr==asym
ans=
1
asym2=asym+sym('0.1')%<12>
astr2=asym+'0.1'%<13>
asym2=
3^(2/3)+0.6
astr2=
3^(2/3)+0.6
(5)
get(0,'format')%<14>
format%<15>
formatcompact%<16>
disp('字符串''0.1''的ASCII码值数组')
disp(double('0.1'))
ad
ans=
long
字符串'0.1'的ASCII码值数组
484649
ad=
2.5801
ad3=ad+double('0.1')%<20>
astr3=ad+'0.1'%<21>
ad3=
50.580148.580151.5801
astr3=
50.580148.580151.5801
(6)
asym+sym(ad)+sym('0.1')%<22>
asym_ad_str=asym+ad+'0.1'%<23>
ans=
3^(2/3)+3.180********19040905619476689026
asym_ad_str=
3^(2/3)+3.1800838230519040905619476689026
(7)
da=double('a')
sda=str2double('a')
d3=double('3')
sd3=str2double('3')
da=
97
sda=
NaN
d3=
51
sd3=
3
103符号数字向双精度数字转换
【例5.2-3】double指令的不同作用。
(1)
clear
formatlong
ad=1/2+sqrt
(2)
astr='1/2+sqrt
(2)'
asym=sym(astr)
ad=
1.914213562373095
astr=
1/2+sqrt
(2)
asym=
2^(1/2)+1/2
(2)
double_asym=double(asym)
double_asym==ad