1、Matlab程序設計實踐與提高（數學研-13，2013秋）第4講符號計算內容簡介：（1）符號對象和符號表達式的創建、操作；（2）微積分、微分方程、矩陣、代數方程的符號計算；（3）符號計算結果的可視化。符號計算工具箱 （Symbolic Math Toolbox ）符號計算引擎：從Matlab 7.9 (2008b) 開始采用 MuPAD（Multi Processing Algebra Data Tool），之前版本采用Maple。符號計算：解算數學表達式、方程時，不是在離散化的數值點上進行，而是憑借一系列恒等式和數學定理，通過推理和演繹，獲得解析結果。這種符號計算建立在數值完全準確表達和推演

2、嚴格解析的基礎之上，因此所得結果是完全準確的。計算機內部是如何處理、理解、推導數學中的“符號”？延伸閱讀推薦：可以看看 計算機代數以及“數學機械化”J4.1 符號對象和符號表達式4.1.1 符號對象的創建和衍生數學表達式和方程的基本組成：（1）數字、參數、變量；（2）運算符號（）；（3）數學函數（如sin、exp等）。（i） 請記住以下兩個Matlab規則：（1）任何基本符號對象（數字、參數、變量、表達式）都必須借助專門的符號函數指令sym或syms定義；（2）任何包含符號對象的表達式或方程，將繼承符號對象的屬性，即任何包含符號對象的表達式、方程也一定是符號對象。（ii） 符號數字一般有限字長

3、數字計算機在記述數字時，不能始終保證完全準確。而在符號計算中，需參與計算的數字完全準確，因此Matlab在符號計算時采用的數字是所謂的符號（類）數字。定義格式為：sym(Num) 創建一個符號數字Num（在符號表達式中使用）sc=sym(Num) 創建一個符號常數sc，該常數值準確等于Num提醒（1）當且僅當Num為整數（如2）、有理數（如4/123）、預定義常數（如pi）時，才能保證：sym指令創建的符號數字（或數字表達式）是完全精準的。 （2）Num應當杜絕使用“十進制小數（如0.123）”、“科學記述數（如1e-2）”，以確保創建的符號數、符號表達式精準。【例4.1-1】演示準確的及應杜

4、絕的符號對象創建。（1）完全精準的符號數字或表達式的創建【強烈推薦的】cleara=sin(3/10) %雙精度的數值類常數sa=sym(sin(3/10) %利用有理數創建一個完全精準的符號數值表達式vpa(sa - a, 40) %計算兩者間的40位精度誤差 a = 0.2955sa =sin(3/10)ans =0.00000000000000002892172294171061964864168785515107729858 （2）不完全精準的符號數字或表達式的創建【應杜絕的】sa1=sym(sin(0.3)sa2=sym(sin(3e-1)vpa(sa-sa2, 40) sa1 =0

5、.29552020666133957510532074568503sa2 =0.29552020666133957510532074568503ans =-6.349448172324523418097794487226753809579e-41 （3）可以觀察一下數據的類型whos Name Size Bytes Class Attributes a 1x1 8 double ans 1x1 60 sym sa 1x1 60 sym sa1 1x1 60 sym sa2 1x1 60 sym 說明前面的提醒（針對數字表達形式的限制），是由符號計算引擎MuPAD的特性決定的。在早期版本的Mat

6、lab中，采用的Maple符號計算引擎則無此要求。（iii） 基本符號變量通常，對于表達式，a、b稱為參數，x稱為變量。在Matlab的符號計算中，a、b、c統稱為“基本符號變量”。 而當對表達式求解、繪圖等操作時，假如不做專門設定，那么x總被默認為“待解符號變量”或稱“自由符號變量”，而其他基本符號變量被作為“符號參數”處理。請help sym查看定義基本符號變量的指令格式：sym(para) syms para 定義單個復數域符號變量parasym(para, flag) syms para flag 定義單個flag指定域符號變量parasyms para1 para2 paraNsym

7、s para1 para2 paraN flag說明（1）若沒有flag指定，則默認為“復數域”；（2）flag： real 實數域，positive 正實數域；（3）在syms para1 para2 中，para1 和 para2 之間用空格分隔；（4）sym和syms兩者功能有異同，。（iv） 自由符號變量基本符號變量可分為：自由符號變量、符號參數。 解題通常是圍繞自由符號變量進行，而解得的結果通常是“用符號參數構成的表達式表述自由符號變量”。在解題時，自由符號變量可由軟件“默認地自動認定”，也可以“人為指定”。可借助symvar指令，來了解該“默認規則”：在沒有專門指定變量名的符號計算

8、中，Matlab將按照與小寫字母x的ASCII碼距離自動識別自由符號變量。 （故而在symvar識別，字母優先次序為 x，y，w，z，v ）help symvar 【例4.1-2】用符號計算研究方程的解。（本例演示：符號常數、基本符號變量、自由符號變量的區別；symvar兩種調用格式的不同功能；Matlab中符號計算默認識別自由符號變量的能力；同一方程，指定的待解變量不同，所得解也不同。）(1)syms u v w z a5f=sym(3);eq = sin(f) * u * z2 + v * z + f*w - a5; (2)symvar(eq) ans = a5, u, v, w, z %

9、 按字母順序列出，注意沒有fsymvar(eq, 100) ans = w, z, v, u, a5 %按離x的距離列出symvar(eq, 1) ans =w (3)result1 = solve(eq) % solve(eq, w)result1 =a5/3 - (v*z)/3 - (u*z2*sin(3)/3 % w = (4) result2 = solve(eq, z) result2 = -(v + (v2 + 4*a5*u*sin(3) - 12*u*w*sin(3)(1/2)/(2*u*sin(3) -(v - (v2 + 4*a5*u*sin(3) - 12*u*w*sin(

10、3)(1/2)/(2*u*sin(3) % z = .【例4.1-3】基本符號變量、符號常數、元符號表達式、衍生符號表達式(1)syms a b x X Y %定義基本符號變量k=sym(3); %定義符號常數z=sym(c*sqrt(d) + y*sin(t); %定義“元”符號表達式expr = a * z * X + (b * x2 + k) * Y; % 構造“衍生”符號表達式 (2)symvar(expr) %列出全部基本符號變量 ans = X, Y, a, b, c, d, t, x, y % 注意：k、z 不在此中(3)symvar(expr,3) %列出自認定的最優先的3個自

11、由符號變量 ans = x, y, t 【例4.1-4】symvar確定自由變量是對整個矩陣進行的。演示：符號矩陣的一種創建方法；矩陣在Matlab中是作為一個整體看待的。syms a b t u v x yA=a+b*x, sin(t)+u; x*exp(-t), log(y)+v A = a + b*x, u + sin(t) x*exp(-t), v + log(y) symvar(A) ans = a, b, t, u, v, x, y symvar(A,1) ans =x % 注意結果不是 x, u; x y，也就是說對于矩陣，symvar對自由變量的認定不是對矩陣元素逐個進行，而是

12、對整個矩陣進行認定的。4.1.2 符號計算中的算符與數值計算中的算符完全相同。說明采用了 重載（overload）技術，可做到這點。（會C+編程的同學，可思考一下C+中的重載）（注： 請help size，看“ Overloaded methods: ”）4.1.3 符號計算中的函數指令Matlab提供的是面向對象的軟件環境。對于不同的數據對象（如數值類和符號類），借助重載技術，把具有相同函數計算功能的文件采用同一個函數名加以保存。這樣處理后，雖然不同類型數據的處理方法不同，但從形式上看，用于數值計算的函數與用于符號計算的函數卻沒有什么區別。至于在運算中，調用哪個函數文件，完全由所計算的對象屬

13、性（數值類還是符號類）決定。Matlab用于符號計算的函數很多（見表4.1），大致可分為三個層次： 第一層次：幾乎與所有數值類函數和指令對應的“同名符號類函數和指令”。 第二層次：約50個經典特殊函數（如誤差函數、貝塞爾函數、橢圓積分等）,它們要借助mfun調用，在Matlab數值計算中沒有對應的函數可供調用。 第三層次：數量很大的MuPAD庫函數，它們借助evalin和feval指令調用。表4.1 MATLAB中可供調用的符號計算函數指令類別情況描述與數值計算對應關系基本函數三角函數、雙曲函數及反函數；除atan2外名稱和使用方法相同指數、對數函數（如exp，expm）名稱和使用方法相同復數

14、函數（注意：沒有幅角函數angle）名稱和使用方法相同矩陣分解函數（如eig、svd）名稱和使用方法相同方程求解函數solve不同微積分函數（如diff，int）不完全相同積分變換和反變換函數（如laplace，ilaplace）只有離散Fourier變換繪圖函數（如ezplot，ezsurf）數值繪圖指令更豐富經典特殊函數如誤差函數erf、貝塞爾函數besselj、第一類完全橢圓積分EllipticK等；通過mfunlist可以看到所有經典函數名無對應函數MuPAD庫函數借助evalin和feval指令可調用比mfunlist所列范圍更廣泛的MuPAD庫函數；需要具備MuPAD語言知識無對應

15、函數說明雖然數值計算和符號計算中有許多同名函數，但是在實際使用函數時，還是要十分注意函數對數據類型的要求，否則容易出錯。比如，對數字而言，有雙精度數字和符號數字之分，這兩種數字在顯示時形式上非常相似，但假如把符號類數字輸入到某個只對數值數據適用的函數（如plot），就一定會產生錯誤。4.1.4 符號對象的識別有的函數適用于多種數據對象，有的只對某種適用。在數值計算和符號計算混合適用情況下，由于函數指令和數據對象不適配引起的錯誤很容易發生。為避免這種錯誤，Matlab提供了用于識別數據對象屬性的指令： class(var) 給出變量var的數據類別（如double，sym等）isa(var, O

16、bj) 若變量var是Obj代表的類別，給出1，表示“真”whos 給出所有Matlab內存變量的屬性【例4.1-5】數據對象及其識別指令的使用。(1) cleara=1; b=2; c=3; d=4;Mn = a, b; c, dMc = a, b; c, dMs = sym(Mc) Mn = 1 2 3 4Mc =a, b; c, dMs = a, b c, d (2)size(Mn)size(Mc)size(Ms) ans = 2 2ans = 1 12ans = 2 2 (3)class(Mn)class(Mc)class(Ms) ans =doubleans =charans =sy

17、m (4)isa(Mn, double)isa(Mc, char)isa(Ms, sym) ans = 1ans = 1ans = 1 (5) whos Mn Mc Ms Name Size Bytes Class Attributes Mc 1x12 24 char Mn 2x2 32 double Ms 2x2 60 sym 4.1.5 符號運算機理和變量假設（i） 符號運算的工作機理Matlab的符號計算是由MuPAD引擎在其專有的內存工作空間中執行，而只是把計算結果送回到Matlab的內存空間。（ii） 對符號變量的限定性假設缺省設置：“復數變量”sym(x) sym(x, real)

18、 sym(x, positive)更多限定：用assume和evalin。比如 assume(x 1)（iii） 清除變量和撤銷假設由于符號變量和其假設存放在不同的內存空間，因此刪除符號變量和撤銷關于變量的假設是兩件需要分別處理的事。clear x 清除Matlab內存中的變量xsyms x clear 撤銷MuPAD內存中對變量x的任何假設，恢復為“復數”變量sym(x, clear) 同上assumptions(x) 顯示符號變量x的限定性假設assumptions 顯示MuPAD內存中已帶限定性假設的全部符號變量reset(symengine) 重啟MuPAD引擎，清空MuPAD內存中所

19、有內容。【例4.1-6】syms對變量所作限定性假設的影響。(1) 在默認的復數域求根syms x clear f = x3 + 4.75*x + 2.5;rf=solve(f, x) % 求出 f=0 的全部根assumptions(x) rf = -1/2 (79(1/2)*i)/4 + 1/4 1/4 - (79(1/2)*i)/4ans = empty sym (2) 求實數根syms x real %設 x 是實數rfr = solve(f, x) %求出 f=0 的實根assumptions(x) %獲悉MuPAD內存中關于x的假設rfr =-1/2ans =x in R_ (3)

20、僅在Matlab內存中刪除x，再解新方程x2+x+5=0clear x %僅清除Matlab內存的變量xsyms xg=x2 + x + 5;rg=solve(g,x) %此前 x為實數的假設仍然起作用，所以導致無解Warning: Explicit solution could not be found. In solve at 169 rg = empty sym (4)須在MuPAD解除x為實數的假設再求解sym(x, clear)rg=solve(g,x) ans =xrg = (19(1/2)*i)/2 - 1/2 - (19(1/2)*i)/2 - 1/2 【例4.1-7】。(1)

21、clear allreset(symengine)Da=1.2; Dw=1/3;syms sa sw sx sy szsyms A B positivesyms C real (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)assumptions ans = 0 A, 0 B,

22、 C in R_ (4)clear Awho Your variables are:B C Da Dw ans sa sw sx sy sz assumptions ans = 0 A, 0 B, C in R_ (5)syms B clearwho Your variables are:B C Da Dw ans sa sw sx sy sz assumptions ans = 0 a. LIMIT(F,a) uses symvar(F) as the independent variable. LIMIT(F) uses a = 0 as the limit point. LIMIT(F,

23、x,a,right) or LIMIT(F,x,a,left) specify the direction of a one-sided limit. Examples: syms x a t h; limit(sin(x)/x) returns 1 limit(x-2)/(x2-4),2) returns 1/4 limit(1+2*t/x)(3*x),x,inf) returns exp(6*t) limit(1/x,x,0,right) returns inf limit(1/x,x,0,left) returns -inf limit(sin(x+h)-sin(x)/h,h,0) re

24、turns cos(x) v = (1 + a/x)x, exp(-x); limit(v,x,inf,left) returns exp(a), 0 【例4.3-1】limit示例：和。syms t x ks=sin(k*t)/(k*t);f=(1-1/x)(k*x); Lsk=limit(s, t, 0) % t趨于0時的極限，limit(s), limit(s,0)Ls1=subs(Lsk,k,1) % 給出k=1時極限的值，即 的值Lsk =1Ls1 = 1 Lfk=limit(f,x,Inf)Lf1=vpa(subs(Lfk, k, sym(-1), 48) % 給出48位精度的自然

25、數Lfk =exp(-k)Lf1 =2.7182818284590452353602874713526624977572470937 【例4.3-2】diff示例：，求。syms a t xf=a, t3; t*cos(x), log(x) f = a, t3 t*cos(x), log(x) dfdx = diff(f, x) d2fdx2 = diff(f, t, 2) d2fdxdt = diff(diff(f, x), t) %兩次調用diff求二階混合導數dfdx = 0, 0 -t*sin(x), 1/xd2fdx2 = 0, 6*t 0, 0d2fdxdt = 0, 0 -sin

26、(x), 0 【例4.3-3】，求。分析：可知在x=處存疑外，f(x) 是處處光滑可導的。(1)對于x0，根據導數定義 求右導數clearsyms xsyms d positiveFp=sin(x); % 在x0時，f(x)=sin(x) p: positivedFp=limit(subs(Fp,x,x+d)-Fp)/d,d,0) %求x0區間的導數dFp0=limit(subs(Fp,x,d)-subs(Fp,x,0)/d,d,0) % 求x=0+的右導數 dFp =cos(x)dFp0 =1 (2)對于x0，根據導數定義 求左導數Fn=sin(-x); % n: negativedFn=l

27、imit(Fn-subs(Fn,x,x-d)/d,d,0) dFn0=limit(subs(Fn,x,0)-subs(Fn,x,-d)/d,d,0) %求x=0-的左導數 dFn =-cos(x)dFn0 =-1 (3) 直接用 diff 求導數f=sin(abs(x);dfdx=diff(f, x)dfdx0=subs(dfdx,x,0) % 根據f(x), x=0處有但dfdx =cos(abs(x)*sign(x)dfdx0 = 0 (4)圖形觀察clf;xn=-3/2*pi:pi/50:0; xp=0:pi/50:3/2*pi; xnp=xn,xp(2:end);figure; hol

28、d on;plot(xnp,subs(f, x, xnp), k, LineWidth, 4) ;plot(xn, subs(dFn, x, xn), -r,LineWidth, 2);plot(xp, subs(dFp, x, xp), -b, LineWidth, 2);legend(char(f), char(dFn), char(dFp);grid on; hold off; 【例4.3-4】設，求。本例演示：如何實現隱函數求導。(1)對方程（即隱函數）求導clearsyms xg=sym(cos(x+sin(y(x)=sin(y(x) % 把y寫成y(x)，表明y是x的函數 dgdx

29、=diff(g,x) g =cos(x + sin(y(x) = sin(y(x)dgdx =-sin(x + sin(y(x)*(cos(y(x)*diff(y(x), x) + 1) = cos(y(x)*diff(y(x), x) (2)用符合規則的新變量名dydx替代dgdx中的 diff(y(x), x)，因后者不是變量名dgdx1 = subs(dgdx, diff(y(x),x), dydx) dgdx1 =-sin(x + sin(y(x)*(dydx*cos(y(x) + 1) = dydx*cos(y(x) (3)對變量dgdx1代表的符號方程，求解dydx，使通過x、y表

30、達出來。dydx = solve(dgdx1, dydx) dydx =-sin(x + sin(y(x)/(cos(y(x) + cos(y(x)*sin(x + sin(y(x) 【例4.3-5】taylor示例：求在處展開的8階Maclaurin級數。(1) 用符號數學工具箱里提供的taylor函數syms xr=taylor(x*exp(x), x, 0, Order, 9)pretty(r) % 展開多項式的易讀格式r =x8/5040 + x7/720 + x6/120 + x5/24 + x4/6 + x3/2 + x2 + x 8 7 6 5 4 3 x x x x x x 2

31、 - + - + - + - + - + - + x + x 5040 720 120 24 6 2 (2)直接調用MuPAD引擎計算R=evalin(symengine, series(x*exp(x), x=0,8)pretty(R) R =Series:Puiseux:create(1, 1, 9, 1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, x, 0, Undirected) 3 4 5 6 7 8 2 x x x x x x 9 x + x + - + - + - + - + - + - + O(x ) 2 6 24 120 720 504

32、0 4.3.2 序列/級數的符號求和、求積s=symsum(f,v,a,b) 求通式f在指定變量v取遍a,b中所有整數時的和，即求help symsum - help for sym/symsum - SYMSUM Symbolic summation. SYMSUM(f) evaluates the sum of a series, where expression f defines the terms of a series, with respect to the default symbolic variable defaultVar determined by symvar. The value of the default variable changes from 0 to defaultVar - 1. SYMSUM(f,x) evaluates the sum of a series, where expression f defines the terms of a series, with respect to the symbolic vari