理学符号计算1.docx

1、理学符号计算1第 5 章 符号计算 符号计算的优点：凭借恒等式，数学定理，通过推理和演绎，给出具有“无限尺度”描写能力的解析结果。当没有封闭解时，符号计算则妥协地给出任意精度数值解。所谓“无限尺度”是相对数值计算的“有限精度、有限空间”而言的。 MATLAB本身由数值计算引擎驱动，而没有符号演绎能力。2008年前，MATLAB的符号计算能力借助于Maple。现在，随MATLAB默认安装的符号计算引擎是MuPAD。虽然MATLAB尽力地保持着形式上的向前兼容，但引擎换装确实导致MATLAB符号计算环境发生了根本改变。本章内容完全针对MuPAD引擎的符号计算而写，分三个层面： 无需任何MuPAD知

2、识，仅使用MATLAB符号数学工具包（Symbolic math toolbox）提供的（前台）函数实施的符号计算及仿真。本章第1节比较完整地描述了符号计算的机理、规则和帮助系统。第2到第6节所涉内容包括微积分、微分方程、积分变换、矩阵分析和代数方程。第7节代数状态方程和第8节数据探索的内容，用以表现现代计算能力对传统方法或技巧的冲击。 需要少量MuPAD知识，拓展符号计算应用范围的内容安排在第5.9节。在这一节中不仅讲述特殊函数计算，而且还介绍了如何利用符号函数产生M函数文件，如何利用符号函数制作用户所需的Simulink模块。 需要较多MuPAD知识，借助evalin和feval指令进入M

3、uPAD空间完成符号计算。 对于那些在MATLAB中仍选择Maple符号计算引擎的读者来说，本章第一层面的内容仍可借鉴。 值得指出：随书光盘mbook目录上的“ch05_符号计算.doc”保存有本章全部算例的运作指令、计算结果和图形；mfile目录上则保存有本章算例中所有带exm前缀文件名的M文件和MEX文件。.1 符号对象的产生和识别.1.1 基本符号对象的创建10 1 定义符号数字和符号常数10 2 定义基本符号变量10 3 定义元符号表达式.1.2 符号计算中的算符和函数指令10 1 符号计算中的算符10 2 符号计算中的函数指令.1.3 符号对象、变量、自由变量的识别10 1 符号对象

4、的识别【例5.1-1】数据对象及其识别指令的使用。（1）cleara=1;b=2;c=3;d=4; Mn=a,b;c,d Mc=a,b;c,d Ms=sym(Mc) Mn = 1 2 3 4Mc =a,b;c,dMs = a, b c, d （2）SizeMn=size(Mn) SizeMc=size(Mc) SizeMs=size(Ms) SizeMn = 2 2SizeMc = 1 9SizeMs = 2 2 （3）CMn=class(Mn)CMc=class(Mc)CMs=class(Ms) CMn =doubleCMc =charCMs =sym （4）isa(Mn,double)is

5、a(Mc,char)isa(Ms,sym) ans = 1ans = 1ans = 1 （5）whos Mn Mc Ms Name Size Bytes Class Attributes Mc 1x9 18 char Mn 2x2 32 double Ms 2x2 60 sym 10 2 符号变量及自由变量的认定【例5.1-2】基本符号变量和衍生符号表达式的定义、符号变量、自由符号变量的机器辨认。（1）clear syms a b k t w x y z X A c=sym(22); f=sym(M*N); ex1=c+f+(exp(-a*t)*sin(w*t)+b*y+k*x+A*X+log

6、(z) ex1 =log(z) + b*y + k*x + sin(t*w)/exp(a*t) + A*X + M*N + 22 （2）symvar(ex1) % 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) % 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 % % Your variables a

7、re:A X a ans b c ex1 f k t w x y z （5）dex1dx=diff(ex1,z) %iex1db=int(ex1,b) % dex1dx =1/ziex1db =(y*b2)/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

8、 + sin(t*w)/exp(a*t) + A*X + M*N + 22 （2）symvar(ex2) % ans = A, M, N, X, a, b, k, t, w, x, y, z （3）symvar(ex2,20)% ans = x, y, w, z, t, k, b, a, X, N, M, A （4）who % Your variables are:ans c ex2 f （5）dex1dx=diff(ex2,z) %? Undefined function or variable z. dex1db=int(ex2,b) % ? Undefined function or v

9、ariable b. （6）E3=sym(a*sqrt(theta) % ? Error using = sym.symconvertExpression at 2547Error: argument must be of Type:Arithmetical sqrtError in = sym.symconvertChar at 2458 s = convertExpression(x);Error in = sym.symconvertCharWithOption at 2441 s = convertChar(x);Error in = sym.symtomupad at 2195 S

10、= convertCharWithOption(x,a);Error in = sym.symsym.sym at 111 S.s = tomupad(x,); E4=sym(a*sqrt(theta123) E4 =a*theta123(1/2) 【例5.1-4】symvar确定自由变量是对整个数组进行的。（1）clearsyms a b t v w x y zA=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

11、) ans = x, y, w symvar(A,10) ans = x, y, w, z, v, t, b, a .1.4 符号运算机理和变量假设10 1 符号运算的工作机理10 2 对符号变量的限定性假设（1）（2）10 3 清除变量和撤销假设【例5.1-5】syms对变量所作限定性假设的影响。（1）syms x clearf=x3+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) % ans =C_ evalin(symengine,

12、property:showprops(x);) % ans = empty sym （2）syms x real rfr=solve(f,x) rfr =-1/2 evalin(symengine,getprop(x) % ans =R_ evalin(symengine,property:showprops(x);)% ans =x in R_ （3）clear x syms xg=x2+x+5;rg=solve(g,x) Warning: Explicit solution could not be found. In solve at 81rg = empty sym （4）syms x

13、 clear rg=solve(g,x) rg = - 1/2 - (19(1/2)*i)/2 - 1/2 + (19(1/2)*i)/2 【例5.1-6】借助evalin对符号变量进行假设的设定和撤销。（1）clearsyms xf=x3-5.25*x-2.5;rf=feval(symengine,solve,f,x) % rf = -1/2, 5/2, -2 （2）x=sym(x,positive); rf=feval(symengine,solve,f,x) % rf =5/2 rfs=solve(f,x) rfs =5/2 evalin(symengine,getprop(x) % a

14、ns =Dom:Interval(0, Inf) evalin(symengine,delete x) evalin(symengine,getprop(x) % ans =C_ （3）evalin(symengine,assume(x,Type:Integer);) %rf=feval(symengine,solve,f,x) % rf =-2 evalin(symengine,assumeAlso(x0) %rf_ev=evalin(symengine,solve(x3-5.25*x-2.5);) %rf_evc=evalin(symengine,solve(,char(f),) %rf_

15、f=feval(symengine,solve,f,x) %rf_ev = empty sym rf_evc = empty sym rf_f = empty sym evalin(symengine,assume(x, Type:Complex);) %rf=feval(symengine,solve,f,x) % rf = -1/2, 5/2, -2 （4）syms % ans f rf rf_ev rf_evc rf_f rfs x .1.5 符号帮助及其他常用指令10 1 符号运作的帮助体系（1）“直接调用符号计算指令”的求助（2）借助mfun调用的“特殊函数指令”的求助（3）借助ev

16、alin或feval调用的“MuPAD指令”的求助。【例5.1-7】关于laplace, erfc, rec三个指令的求助过程。（1）图5.1-1 “可直接调用的符号计算指令”的帮助信息搜索（2）图5.1-2 MATLAB帮助浏览器展示全部特殊函数符号计算指令及用法（3）图5.1-3 MuPAD帮助浏览器展示的rec指令帮助信息10 2 服务于符号运作的其他指令【例5.1-8】各种帮助指令、符号变量罗列指令（1）clear all reset(symengine) %Da=1.2;Dw=1/3; %syms sa sw sx sy sz %syms A B positive %syms C r

17、eal % （2）whos % Name Size Bytes Class Attributes A 1x1 60 sym B 1x1 60 sym C 1x1 60 sym Da 1x1 8 double Dw 1x1 8 double sa 1x1 60 sym sw 1x1 60 sym sx 1x1 60 sym sy 1x1 60 sym sz 1x1 60 sym （3）syms % A B C sa sw sx sy sz （4）evalin(symengine,anames(Properties) % ans = A, B, C （5）clear A syms % B C an

18、s sa sw sx sy sz evalin(symengine,anames(Properties) % ans = A, B, C （6）syms B clear %syms B C ans sa sw sx sy sz evalin(symengine,anames(Properties) % ans = A, C .2 数字类型转换及符号表达式操作.2.1 数字类型及转换10 1 三种数字类型及转换指令图 5.2-1 符号、字符、数值间的相互转换10 2 双精度数字向符号数字转换【例5.2-1】本例目的：准确符号数字的产生；双精度数字转换成符号数字的各种格式；由双精度数字转换而得符号

19、数的误差；vpa的用法；数字类别的判断。（1）clear sa=sym(1/3) %sb=sym(pi+sqrt(5) % sa =1/3sb =pi + 5(1/2) （2）format long a=1/3 %b=pi+sqrt(5) %sa_a=vpa(sa-a) %sb_b=vpa(sb-b) %a = 0.333333333333333b = 5.377660631089583sa_a =0.0sb_b =0.000000000000000013822375841085200048593542564188 （3）asr=sym(a) bsr=sym(b) sa_asr=vpa(sa-

20、asr)sb_bsr=vpa(sb-bsr) asr =1/3bsr =189*1719/35184372088832sa_asr =0.0sb_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/12bse =189*1719/35184372088832sa_ase =0.083333333333333333333333333333333*epssb_bse =0.000

21、000000000000013822375841085200048593542564188 （5）asf=sym(a,f) bsf=sym(b,f) sa_asf=vpa(sa-asf)sb_bsf=vpa(sb-bsf) asf =6004799503160661/180*481984bsf =189*1719/35184372088832sa_asf =0.000000000000000018503717077085942340393861134847sb_bsf =0.000000000000000013822375841085200048593542564188 （6）asd=sym(

22、a,d) bsd=sym(b,d) sa_asd=vpa(sa-asd)sb_bsd=vpa(sb-bsd) asd =0.33333333333333331482961625624739bsd =5.3776606310895829210494412109256sa_asd =0.000000000000000018503717077085943333333327037792sb_bsd =0.000000000000000013822375841085179119638453173161 （7）class(sa) ans =sym disp(a的类别 asr的类别 sa_a的类别) dis

23、p(class(a),blanks(4),class(asr),blanks(7),class(sa_a) a的类别 asr的类别 sa_a的类别double sym sym 【例5.2-2】本例目的：揭露字符串数字在不同计算引擎中的不同表现；指令format的用法；double与str2double的区别。（1）clearformat long %ad=1/2+3(2/3) astr=1/2+3(2/3) asym=sym(astr) ad = 2.580083823051904astr =1/2+3(2/3)asym =3(2/3) + 1/2 （2）disp( ad astr asym)

24、disp(blanks(6), 的数据类别) disp(class(ad),blanks(3),class(astr),blanks(4),class(asym) ad astr asym 的数据类别double char sym （3）asymPLUSad=asym+sym(0.1) %asymPLUSda2=asym+0.1 % asymPLUSad =3(2/3) + 3/5asymPLUSda2 =3(2/3) + 3/5 （4）astr=asym ans = 1 asym2=asym+sym(0.1) %astr2=asym+0.1 %asym2 =3(2/3) + 0.6astr2

25、 =3(2/3) + 0.6 （5）get(0,format) %format %format compact %disp(字符串0.1的ASCII码值数组)disp(double(0.1) ad ans =long字符串0.1的ASCII码值数组 48 46 49ad = 2.5801 ad3=ad+double(0.1) %astr3=ad+0.1 % ad3 = 50.5801 48.5801 51.5801astr3 = 50.5801 48.5801 51.5801 （6）asym+sym(ad)+sym(0.1) %asym_ad_str=asym+ad+0.1 %ans =3(2

26、/3) + 3.180*19040905619476689026asym_ad_str =3(2/3) + 3.1800838230519040905619476689026 （7）da=double(a) sda=str2double(a) d3=double(3) sd3=str2double(3) da = 97sda = NaNd3 = 51sd3 = 3 10 3 符号数字向双精度数字转换【例5.2-3】double指令的不同作用。（1）clearformat long ad=1/2+sqrt(2) astr=1/2+sqrt(2) asym=sym(astr) ad = 1.914213562373095astr =1/2+sqrt(2)asym =2(1/2) + 1/2 （2）double_asym=double(asym) double_asym=ad double_asym =