理学符号计算1.docx
《理学符号计算1.docx》由会员分享，可在线阅读，更多相关《理学符号计算1.docx（41页珍藏版）》请在冰豆网上搜索。
理学符号计算1
第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
double_asym=
升级会员 