1、MatLab & 數學建模授課:唐 靜 波 (九江學院理學院第四講 數值計算符號數學工具箱符號表達式的運算numeric 符號到數值的轉換pretty 顯示悅目的符號輸出subs 替代子表達式sym 建立符號矩陣或表達式symadd 符號加法symdiv 符號除法symmul 符號乘法symop 符號運算sympow 符號表達式的冪運算symrat 有理近似symsub 符號減法symvar 求符號變量符號表達式的簡化collect 合并同類項expand 展開factor 因式simple 求解最簡形式simplify 簡化symsum 和級數符號多項式charpoly 特征多項式h

2、orner 嵌套多項式表示numden 分子或分母的提取poly2sym 多項式向量到符號的轉換sym2poly 符號到多項式向量的轉換符號微積分diff 微分int 積分jordan 約當標準形taylor 泰勒級數展開符號可變精度算術digits 設置可變精度vpa 可變精度計算求解符號方程compose 函數的復合dsolve 微分方程的求解finverse 函數逆linsolve 齊次線性方程組的求解solve 代數方程的求解符號線性代數charploy 特征多項式determ 矩陣行列式的值eigensys 特征值和特征向量inverse 矩陣逆jordan 約當標準形linsolv

3、e 齊次線性方程組的解transpose 矩陣的轉置一、 方程求解求解單個代數方程MATLAB 具有求解符號表達式的工具 , 如果表達式不是一個方程式 (不含等 號 ,則在求解之前函數 solve 將表達式置成等于 0。>> solve( ' a*x2+b*x+c ' % solve for the roots of the eqution ans=1/2/a*(-b+(b2-4*a*c1/21/2/a*(-b-(b2-4*a*c1/2結果是符號向量, 其元素是方程的 2個解。 如果想對非缺省 x 變量求解, solve 必須指定變量。>> solve(

4、 ' a*x2+b*x+c ' , ' b ' % solve for b ans=-(a*x2+c/x帶有等號的符號方程也可以求解。>> f=solve( ' cos(x=sin(x ' % solve for xf=1/4*pi>> t=solve( ' tan(2*x=sin(x ' t= 0acos(1/2+1/2*3(1/2acos(1/2=1/2*3(1/2并得到數值解。>> numeric(fans=0.7854>> numeric(tans=0 + 0.8314i1.

5、9455注意在求解周期函數方程時,有無窮多的解。在這種情況下, solve 對解的 搜索范圍限制在接近于零的有限范圍,并返回非唯一的解的子集。如果不能求得符號解,就計算可變精度解。>> x=solve( ' exp(x=tan(x ' x=1.306326940423079代數方程組求解可以同時求解若干代數方程,語句 solve (s1, s2, . , sn 對缺省變量求 解 n 個方程,語句 solve(s1, s2, . , sn , ' v1, v2, . , vn '對 n 個 ' v1, v2, .vn ' 的未知數求解

6、n 個方程。solve (f 解符號方程式 f 。solve(f1,fn 解由 f1,fn 組成的聯立方程式。我們先定義以下的方程式:>>eq1 = 'x-3=4' % 注意也可寫成 'eq1=x-7'>>eq2 = 'x*2-x-6=0' % 注意也可寫成 'eq2=x*2-x-6'>>eq3 = 'x2+2*x+4=0'>>eq4 = '3*x+2*y-z=10'>>eq5 = '-x+3*y+2*z=5'>>

7、;eq6 = 'x-y-z=-1'>>solve(eq1ans=7>>solve(eq2ans=3,-2' % 原方程式有二個根 3, -2>>solve(eq3ans=-1+i*3(1/2,-1-i*3(1/2' % 注意實根和虛根的表示式>>solve(eq4,eq5,eq6 % 解三個聯立方程式ans=x = -2, y = 5, z = -6如何處理中小學典型的代數問題?黛安娜 (Diane想去看電影,她從小豬存錢罐倒出硬幣并清點,她發現:10美分的硬幣數加上 5美分的硬幣總數的一半等于 25美分的硬幣數。

8、 1美分的硬幣數比 5美分、 10美分以及 25美分的硬幣總數多 10。25美分和 10美分的硬幣總數等于 1美分的硬幣數加上 1/4的 5美分的硬幣 數25美分的硬幣數和 1美分的硬幣數比 5美分的硬幣數加上 8倍的 10美分的 硬幣數多 1。如果電影票價為 3.00美元,爆米花為 1.00美元,糖棒為 50美分,她有足夠的 錢去買這三樣東西?首先,根據以上給出的信息列出一組線性方程,假如 p , n , d 和 q 分別表示 1美分, 5美分, 10美分,和 25美分的硬幣數d n pq p n d q q d p n q p n d +=+-+=+=+-210481然后,建立 MATLA

9、B 符號方程并對變量求解。>> eq1= ' d+(n+p/2=q ' ;>> eq2= ' p=n+d+q-10 ' ;>> eq3= ' q+d=p+n/4 ' ;>> eq4= ' q+p=n+8*d-1 ' ;>>pennies, nickles , dimes , quarters=solve(equ1, equ2, equ3, equ4, ' p, n , d , q ' pennies=16nickles=8dimes=3quarters=1

10、5所以,黛安娜有 16枚 1美分的硬幣, 8枚 5美分的硬幣, 3枚 10美分的硬幣, 15枚 25美分的硬幣,這就意味著>> money=.01*16+.05*8+.10*3+.25*15money=4.6100她就有足夠的錢去買電影票,爆米花和糖棒并剩余 11美分。【例】求解二元函數方程組 =+=-=0 cos( , (0 sin( , (21y x y x f y x y x f 的零點。 (0從三維坐標初步觀察兩函數圖形相交情況x=-2:0.05:2;y=x;X,Y=meshgrid(x,y; %產生 x-y 平面上網點坐標 F1=sin(X-Y;F2=cos(X+Y;F0

11、=zeros(size(X;surf(X,Y,F1,xlabel('x',ylabel('y',view(-31,62,hold on,surf(X,Y,F2,surf(X,Y,F0,shading interp,hold off 圖 5.6.3-0 兩函數的三維相交圖(1在某區域觀察兩函數 0等位線的交點情況clear;x=-2:0.5:2;y=x;X,Y=meshgrid(x,y; %產生 x-y 平面上網點坐標F1=sin(X-Y;F2=cos(X+Y;v=-0.2, 0, 0.2; %指定三個等位值,是為了更可靠地判斷 0等位線的存在。 contour(

12、X,Y,F1,v %畫 F1的三條等位線。hold on,contour(X,Y,F2,v,hold off %畫 F2的三條等位線。 圖 5.6.3-1 兩個二元函數 0等位線的交點圖(2從圖形獲取零點的初始近似值在圖 5.6.3-1中,用 ginput 獲取兩個函數 0等位線(即三線組中間那條線 交點的坐標。x0,y0=ginput(2; %在圖上取兩個點的坐標disp(x0,y0-0.7926 -0.78430.7926 0.7843(3利用 fsolve 求精確解。以求(0.7926,7843附近的解為例。本例直接用字符串表達被解函數。 注意:在此, 自變量必須寫成 x(1, x(2。

13、 假如寫成 xy(1, xy(2,指令運行將出錯。fun='sin(x(1-x(2,cos(x(1+x(2' %<12> xy=fsolve(fun,x0(2,y0(2%<13>xy =0.7854 0.7854(4檢驗fxy1=sin(xy(1-xy(2;fxy2=cos(xy(1+xy(2;disp(fxy1,fxy21.0e-006 *-0.0994 0.2019說明指令 <12><13>可用以下任何一組指令取代。(A 內聯函數形式指令fun=inline('sin(x(1-x(2, cos(x(1+x(2'

14、, 'x' %項 'x' 必須有。 xy=fsolve(fun,x0(2, y0(2;(B M 函數文件 形式及指令先用如下 fun.m 表示被解函數(并在搜索路徑上fun.mfunction ff=fun(xff(1=sin(x(1-x(2;ff(2=cos(x(1+x(2;然后運行指令 xy=fsolve('fun',x0(2,y0(2 。第四步檢驗中的結果表明:所找零點處的函數值小于 610 ,是一個十分接近 零的小數。該精度由 options.TolFun 控制。 options.TolFun 的缺省值是 1.0000e-006。它可以用

15、下列指令看到options=optimset('fsolve'options.TolFunans =1.0000e-006線性方程求解a= 7 2 1 -29 15 3 -2-2 -2 11 51 3 2 13b=4 7 -1 0'x=abx =0.49790.14450.0629-0.0813單個微分方程常微分方程有時很難求解, MATLAB 提供了功能強大的工具,可以幫助求解微 分方程。 函數 dsovle 計算常微分方程的符號解。 因為我們要求解微分方程, 就需 要用一種方法將微分包含在表達式中。 所以, dsovle 句法與大多數其它函數有一 些不同,用字母 D

16、 來表示求微分, D2, D3等等表示重復求微分,并以此來設定方 程。任何 D 后所跟的字母為因變量。MATLAB 解常微分方程式的語法是 dsolve('equation','condition',其中 equation 代表常微分方程式即 y' =g(x , y ,且須 以 Dy 代表一階微分項 y' D2y 代 表二階微分項 y'' , condition 則為初始條件。方程 d y dx22/=0用符號表達式 D2y=0來表示。獨立變量可以指定或由 symvar 規則選定為缺省。例如,一階方程 dy/dx=1+y2的通解為

17、:>> dsolve( ' Dy=1+y2 ' % find the general solution ans=-tan(-x+C1其中, C1是積分常數。求解初值 y(0=1的同一個方程就可產生:>> dsolve(' Dy=1+y2 ', ' y(0=1 ' % add an initial conditiony=tan(x+1/4*pi獨立變量可用如下形式指定:>> dsolve(' Dy=1+y2 ' , ' y(0=1 ' , ' v ' % find

18、solution to dy/dvans=tan(v+1/4*pi讓我們舉一個二階微分方程的例子,該方程有兩個初始條件: d ydx 22=cos(2x-y dy dx(0=0 y(0=1>> y=dsolve(' D2y=cos(2*x-y ', ' Dy(0=0 ', ' y(0=1 ' y=-2/3*cos(x2+1/3+4/3*cos(x>> y=simple(y % y looks like it can be simplified y=-1/3*cos(2*x+4/3*cos(x通常,要求解的微分方程含有一階以

19、上的項,并以下述的形式表示:d ydx 22-2dy dx -3y=0通解為:>> y=solve( 'D2y-2Dy-3*y=0 'y=C1*exp(-x+C2*exp(3*x加上初始條件:y(0=0和 y(1=1可得到:>> y=solve( ' D2y-2Dy-3*y=0 ' , ' y(0=0, y(1=1 ' y=1/(exp(-1-exp(3*exp(-x-1/(exp(-1-exp(3*exp(3*x>> y=simple(y % this looks like a candidate for s

20、implificationy=-(exp(-x-exp(3*x/(exp(3-exp(-1>> pretty(y % pretty it upexp(-x-exp(3 x- -exp(3 -exp(-1現在來繪制感興趣的區域內的結果。>> ezplot(y, -6 2例 :假設有以下三個一階常微分方程式和其初始條件 y' =3x 2, y (2=0.5y' =2. x . cos(y2, y (0=0.25y' =3y+exp(2x, y(0=3對應上述常微分方程式的符號運算式為:>>soln_1 = dsolve('Dy =

21、 3*x2','y(2=0.5' ans=x3-7.500000000000000>>ezplot(soln_1,2,4 % 看看這個函數的長相 >>soln_2 = dsolve('Dy = 2*x*cos(y2','y(0 = pi/4'ans=atan(x2+1>>soln_3 = dsolve('Dy = 3*y + exp(2*x',' y(0 = 3'ans=-exp(2*x+4*exp(3*x微分方程組函數 dsolve 也可同時處理若干個微分方程式,下面有

22、兩個線性一階方程。 dx df=3f+4g dg dx =-4f+3g通解為:>> f, g=dsolve( ' Df=3*f+4*g ' , ' Dg=-4*f+3*g ' f=C1*exp(3*x*sin(4*x+C2*exp(3*x*cos(4*xg=-C2*exp(3*x*sin(4*x+C1*exp(3*x*cos(4*x加上初始條件:f(0=0和 g(0=1,我們可以得到:>> f, g=dsolve( ' Df=3*f+4*g ' , ' Dg=-4*f+3*g ' , ' f(0=0

23、, g(0=1 ' f=exp(3*x*sin(4*xg=exp(3*x*cos(4*x微分和積分微分和積分是微積分學研究和應用的核心,并廣泛地用在許多工程學科。 MATLAB 符號工具能幫助解決許多這類問題。微分符號表達式的微分以四種形式利用函數 diff :>> f= ' a*x3+x2-b*x-c ' % define a symbolic expressionf=a*x3+x2-b*x-c>> diff(f % differentiate with respect to the default variable xans=3*a*x2+2

24、*x-b>> diff(f, 'a ' % differentiate with respect to a ans=x3>> diff(f, 2 % differentiate twice with respect to x ans=6*a*x+2>> diff(f, ' a ' , 2 % differentiate twice with respect to aans=函數 diff 也可對數組進行運算。 如果 F 是符號向量或數組, diff(F對數組內 的各個元素進行微分。>> F=sym(' a*

25、x, b*x2; c*x3, d*s ' % create a symbolic arrayF= a*x, b*x2c*x3, d*s>> diff(F % differentiate the element with respect to x ans= a, 2*b*x3*c*x2, 0注意函數 diff 也用在 MATLAB ,計算數值向量或矩陣的數值差分。對于一個數 值向量或矩陣 M , diff(M計算 M(2: m, : -M(1: m-1, : 的數值差分,如下 所示:>> m=(1: 8.2 % create a vectorM=1 4 9 16

26、25 36 49 64>> diff(M % find the differences between elements ans=3 5 7 9 11 13 15如果 diff 的表達式或可變參量是數值, MATLAB 就非常巧妙地計算其數值差 分;如果參量是符號字符串或變量, MATLAB 就對其表達式進行微分。積分積分 函數 int(f,其中 f 是一符號表達式,它力圖求出另一符號表達式 F 使 diff(F=f。正如從研究微分學所了解的,積分比微分復雜得多。積分或逆求導 不一定是以封閉形式存在, 或許存在但軟件也許找不到, 或者軟件可明顯地求解, 但超過內存或時間限制。 當

27、MATLAB 不能找到逆導數時, 它將返回未經計算的命令。>> int( ' log(x/exp(x2 ' % attempt to integrate ans=log(x/exp(x2同微分一樣, 積分函數有多種形式。 形式 int(f相對于缺省的獨立變量求逆 導數;形式 (f, ' s '相對于符號變量 s 積分;形式 int(f, a , b 和 int(f, ' s ' , a , b , a , b 是數值,求解符號表達式從 a 到 b 的定積分;形式 int(f, ' m ' , ' n '

28、和形式 int(f, ' s ', ' m ', ' n ',其中 m , n 是符號變量,求解符號 表達式從 m 到 n 的定積分。>> f=' sin(s+2*x ' % crate a symbolic functionf=sin(s+2*x>> int(f % integrate with respect to xans=-1/2*cos(s+2*x>> int(f, ' s ' % integrate with respect to sans=-cos(s+2*x>

29、;> int(f, pi/2, pi % integrate with respect to x from/2 toans=-cos(x>> int(f, ' s ', pi/2, pi % integrate with respect to s from /2 to ans=cos(2*x-sin(2*x>> int(f, ' m ', ' n ' % integrate with respect to x from m to nans=-1/2*cos(s+2*n+1/2*cos(s+2*m正如函數 diff 一

30、樣,積分函數 int 對符號數組的每一個元素進行運算。>> F=sym( ' a*x, b*x2; c*x3, d*s ' % create a symbolic arrayF= a*x, b*x2c*x3, d*s>> diff(F % ubtegrate the array elements with respect to x ans=1/2*a*x2, 1/3*b*x31/4*c*x4, d*s*xdiff 函數用以演算一函數的微分項,相關的函數語法有下列 4個:diff(f 傳回 f 對預設獨立變數的一次微分值diff(f,'t'

31、傳回 f 對獨立變數 t 的一次微分值diff(f,n傳回 f 對預設獨立變數的 n 次微分值diff(f,'t',n傳回 f 對獨立變數 t 的 n 次微分值先定義下列三個方程式,接著再演算其微分項: >>S1 = '6*x3-4*x2+b*x-5'>>S2 = 'sin(a'>>S3 = '(1 - t3/(1 + t4'>>diff(S1ans=18*x2-8*x+b>>diff(S1,2ans=36*x-8>>diff(S1,'b'an

32、s=x>>diff(S2ans=cos(a>>diff(S3ans=-3*t2/(1+t4-4*(1-t3/(1+t42*t3 >>simplify(diff(S3ans=t2*(-3+t4-4*t/(1+t42int 函數用以演算一函數的積分項, 這個函數要找出一符號式 F 使得diff(F=f。如果積分式的解析式 (analytical form, closed form 不存在的 話或是 MATLAB 無法找到, 則 int 傳回原輸入的符號式。 相關的函數語法有下列 4個:int(f傳回 f 對預設獨立變數的積分值int(f,'t'傳

33、回 f 對獨立變數 t 的積分值int(f,a,b傳回 f 對預設獨立變數的積分值,積分區間為 a,b, a 和 b 為數值 式int(f,'t',a,b傳回 f 對獨立變數 t 的積分值,積分區間為 a,b, a 和 b 為數 值式int(f,'m','n'傳回 f 對預設變數的積分值,積分區間為 m,n, m 和 n 為符號 式我們示范幾個例子:>>S1 = '6*x3-4*x2+b*x-5'>>S2 = 'sin(a'>>S3 = 'sqrt(x'>&

34、gt;int(S1ans=3/2*x4-4/3*x3+1/2*b*x2-5*x>>int(S2ans=-cos(a>>int(S3ans=2/3*x(3/2>>int(S3,'a','b'ans=2/3*b(3/2- 2/3*a(3/2>>int(S3,0.5,0.6ans=2/25*15(1/2-1/6*2(1/2>>numeric(int(S3,0.5,0.6 % 使用 numeric 函數可以計算積分的數值ans=0.0741數值積分先考慮一個積分式的數學式如下: 其中 a, b分別為這個積分式的

35、上限及下限, f (x 則為要積分的函樹。 要求解上述 的積分式,必須設定 a, b和 f (x 。以 MATLAB 的積分函數求解的過程亦同, 也要定義 f (x 及設定 a,b , 還須設定在區間 a,b 之間離散 點 (discretized points 數目,剩下的工作就是選擇精度不同的積分法來求解了。梯形法MATLAB 提供最簡單的積分函數是梯形法 trapz ,我們先說明梯形法語法 trapz(x,y,其中 x,y 分別代表數目相同 的陣列或矩陣 ,而 y 與 x 的關系可以由 是一函數型態(如 y =sin(x 或是不以函數描述的離散型態。我們看一簡單積分式 以下為 MATLAB 的程式>> x=0:pi/100:pi;>&