MatLab4數(shù)值計(jì)算一

1、MatLab數(shù)學(xué)建模22授課：曾景峰（江西環(huán)境工程職業(yè)學(xué)院）第四講數(shù)值計(jì)算符號數(shù)學(xué)工具箱 符號表達(dá)式的運(yùn)算nu meric符號到數(shù)值的轉(zhuǎn)換p retty顯示悅目的符號輸出subs替代子表達(dá)式sym建立符號矩陣或表達(dá)式symadd符號加法symdiv符號除法symmul符號乘法symop符號運(yùn)算sympow符號表達(dá)式的幕運(yùn)算symrat有理近似symsub符號減法symvar求符號變量符號表達(dá)式的簡化collect expand factor simp le sim plify symsum合并同類項(xiàng)展開因式求解最簡形式 簡化和級數(shù)符號多項(xiàng)式char poly hornernu mde npo

2、Iy2sym sym2 poly特征多項(xiàng)式嵌套多項(xiàng)式表示分子或分母的提取多項(xiàng)式向量到符號的轉(zhuǎn)換 符號到多項(xiàng)式向量的轉(zhuǎn)換符號微積分diffint微分積分jorda n taylor約當(dāng)標(biāo)準(zhǔn)形 泰勒級數(shù)展開符號可變精度算術(shù)digits設(shè)置可變精度vpa可變精度計(jì)算求解符號方程compose函數(shù)的復(fù)合dsolve微分方程的求解fin verse函數(shù)逆lin solve齊次線性方程組的求解Solve代數(shù)方程的求解符號線性代數(shù)char ploy特征多項(xiàng)式determ矩陣行列式的值eige nsys特征值和特征向量in verse矩陣逆jorda n約當(dāng)標(biāo)準(zhǔn)形lin solve齊次線性方程組的解trans

3、pose矩陣的轉(zhuǎn)置方程求解求解單個(gè)代數(shù)方程MATLA具有求解符號表達(dá)式的工具,如果表達(dá)式不是一個(gè)方程式（不含等 號），則在求解之前函數(shù)solve將表達(dá)式置成等于0。 solve( a*xA2+b*x+c )% solve for the roots of the equti onans=1/2/a*(-b+(bA2-4*a*c)A1/2)1/2/a*(-b-(bA2-4*a*c)A1/2)結(jié)果是符號向量，其元素是方程的2個(gè)解。如果想對非缺省x變量求解，solve 必須指定變量。 solve( a*x2+b*x+c ,b ) % solve for bans=-(a*x2+c)/x帶有等號的符號

4、方程也可以求解。 f=solve( cos(x)=sin(x) ) % solve for x f=1/4*pi t=solve( tan(2*x)=sin(x) )t= 0acos(1/2+1/2*3A(1/2)acos(1/2=1/2*3A(1/2)并得到數(shù)值解。 numeric (f) ans=0.7854 numeric(t) ans=00 + 0.8314i1.9455注意在求解周期函數(shù)方程時(shí)，有無窮多的解。在這種情況下， solve 對解的 搜索范圍限制在接近于零的有限范圍，并返回非唯一的解的子集。如果不能求得符號解，就計(jì)算可變精度解。 x=solve( exp(x)=tan(x)

5、 ) x=1.306326940423079代數(shù)方程組求解 的未知數(shù)求解 n 個(gè)方程?？梢酝瑫r(shí)求解若干代數(shù)方程，語句 solve （s1 ，s2， ，sn） 對缺省變量求 解n個(gè)方程，語句 solve（s1 , s2, . , sn, v1 , v2, . , vn ）對n個(gè)v1 , v2， .vnsolve (f)解符號方程式 f 。solve(f1,fn）解由f1,fn組成的聯(lián)立方程式。我們先定義以下的方程式：eq1 = x-3=4; % 注意也可寫成 eq1=x-7eq2 = x2-x-6=0; %注意也可寫成eq2=x*2-x-6eq3 = x2+2*x+4=0;eq4 = 3*x+2

6、*y-z=10;eq5 = -x+3*y+2*z=5;eq6 = x-y-z=-1;solve(eq1) ans=solve(eq2)ans=3,-2 %原方程式有二個(gè)根3, -2solve(eq3) ans=注意實(shí)根和虛根的表示式-1+i*3A(1/2),-1-i*3A(1/2) % solve(eq4,eq5,eq6) % 解三個(gè)聯(lián)立方程式ans= x = -2, y = 5, z = -6如何處理中小學(xué)典型的代數(shù)問題？黛安娜(Diane)想去看電影，她從小豬存錢罐倒出硬幣并清點(diǎn)，她發(fā)現(xiàn):10美分的硬幣數(shù)加上5美分的硬幣總數(shù)的一半等于25美分的硬幣數(shù)。1美分的硬幣數(shù)比5美分、10美分以及2

7、5美分的硬幣總數(shù)多10。25美分和10美分的硬幣總數(shù)等于1美分的硬幣數(shù)加上1/4的5美分的硬幣*硬幣數(shù)多1。25美分的硬幣數(shù)和1美分的硬幣數(shù)比5美分的硬幣數(shù)加上8倍的10美分的如果電影票價(jià)為3.00美元，爆米花為1.00美元，糖棒為50美分，她有足夠的 錢去買這三樣?xùn)|西？首先，根據(jù)以上給出的信息列出一組線性方程, 美分，5美分，10美分，和25美分的硬幣數(shù)假如p, n, d和q分別表示1n + pd+=q p = n +d+q10 q+d2q + p = n +8d 1然后，建立MATLAB號方程并對變量求解。 eq1= d+(n+p)/2=q eq2= p=n+d+q-10 eq3= q+d

8、 二p+n/4 eq4= q+p=n+8*d-1 pennies , nickles , dimes, quarters=solve(equ1 p , n, d, q )penni es=16ni ckles=8dimes=3quarters=15,equ2，equ3, equ4,所以，黛安娜有16枚1美分的硬幣，8枚 5美分的硬幣, 枚25美分的硬幣，這就意味著3枚10美分的硬幣，15 mon ey=01*16+.05*8+.10*3+.25*15mon ey=4.6100她就有足夠的錢去買電影票，爆米花和糖棒并剩余 11美分?！纠壳蠼舛瘮?shù)方程組Jf1（x, ysin（Xy） =0的零

9、點(diǎn)。lf2（x,y） = cos（x + y） =0滬生x-y平面上網(wǎng)點(diǎn)坐標(biāo)(0)從三維坐標(biāo)初步觀察兩函數(shù)圖形相交情況 x=-2:0.05:2;y=x;X, Y=meshgrid(x,y);F1=si n(X-Y );F2=cos(X+Y); F0=zeros(size(X); surf(X ,Y ,F1), xlabel(x),ylabel(y), view(-31,62),hold on, surf(X,Y,F2),surf(X,Y,F0), shadi ng interp, hold off22y-2-20、圖 563-0兩函數(shù)的三維相交圖(1)在某區(qū)域觀察兩函數(shù)0等位線的交點(diǎn)情況cle

10、ar;x=-2:0.5:2;y=x;X,Y=meshgrid(x,y);F1=si n(X-Y);F2=cos(X+Y);滬生x-y平面上網(wǎng)點(diǎn)坐標(biāo)v=-0.2, 0, 0.2;%旨定三個(gè)等位值，是為了更可靠地判斷0等位線的存在。con tour(X,Y,F1,v)hold on,con tour(X,Y,F2,v),hold off%畫 F1的三條等位線。%畫 F2的三條等位線。圖5.6.3-1 兩個(gè)二元函數(shù)0等位線的交點(diǎn)圖(2) 從圖形獲取零點(diǎn)的初始近似值 在圖563-1中，用ginput獲取兩個(gè)函數(shù)0等位線(即三線組中間那條線)交點(diǎn)的坐標(biāo)。%在圖上取兩個(gè)點(diǎn)的坐標(biāo)x0,y0=g inp ut

11、 (2);dis p( x0,y0)-0.7926 -0.78430.79260.7843(3) 利用fsolve求精確解。以求(0.7926,7843 )附近的解為例。本例直接用字符串表達(dá)被解函數(shù)。注意：在此，自變量必須寫成x(1), x(2)。 假如寫成xy(1), xy(2)，指令運(yùn)行將出錯(cuò)。%fun=si n(x(1)-x (2),cos(x(1)+x(2);xy=fsolve(fu n,x0(2),y0 (2)%xy =0.78540.7854(4) 檢驗(yàn)fxy1=si n(xy(1)-xy(2);fxy2=cos(xy(1)+xy (2);dis p(fxy1,fxy2)1.0e-

12、006 *-0.09940.20191說明指令可用以下任何一組指令取代。(A)內(nèi)聯(lián)函數(shù)形式指令%項(xiàng)x必須有。fun=i niin e(s in (x(1)-x (2), cos(x(1)+x (2), x);xy=fsolve(fu n,x0(2), y0(2);(B) M函數(shù)文件形式及指令先用如下fun.m表示被解函數(shù)(并在搜索路徑上)fun .mfunction ff=fun(x)ff(1)=si n(x(1)-x (2);ff(2)=cos(x(1)+x (2);然后運(yùn)行指令 xy=fsolve(fun,x0(2),y0(2)。第四步檢驗(yàn)中的結(jié)果表明：所找零點(diǎn)處的函數(shù)值小于10，是一個(gè)十

13、分接近零的小數(shù)。該精度由 op tio ns.ToIFu n 控制。op tio ns.ToIFu n 的缺省值是 1.0000e-006。它可以用下列指令看到op tio ns=op timset(fsolve);op tio ns.Tol Funans =1.0000e-006線性方程求解a= 7 2 1-215-2-2-21113b=4 7 -1 0x=ab0.49790.14450.0629-0.0813單個(gè)微分方程常微分方程有時(shí)很難求解，MATLAB供了功能強(qiáng)大的工具，可以幫助求解微 分方程。函數(shù)dsovle計(jì)算常微分方程的符號解。因?yàn)槲覀円蠼馕⒎址匠?，就?要用一種方法將微分包含

14、在表達(dá)式中。 些不同，用字母D來表示求微分，D2, 程。任何D后所跟的字母為因變量。MATLA解常微分方程式的語法是所以， dsovle 句法與大多數(shù)其它函數(shù)有一D3等等表示重復(fù)求微分，并以此來設(shè)定方dsolve(equation,condition)equation代表常微分方程式即y =g(x, y)，且須 以Dy代表一階微分項(xiàng)y 表二階微分項(xiàng) y， condition 則為初始條件。，其中D2y代方程d2y/dx2= 0用符號表達(dá)式D2y=(來表示。獨(dú)立變量可以指定或由symvar規(guī)則選定為缺省。例如，一階方程dy/dx=1+y2的通解為: dsolve( Dy=1+yA2 ) % fi

15、nd the gen eral solutionans=-tan (-X+C1)其中，C1是積分常數(shù)。求解初值y(0)=1的同一個(gè)方程就可產(chǎn)生:,y(0)=1 ) % add an initial，y(0)=1 ， v ) % find dsolveC Dy=1+yA2 con diti ony=ta n(x+1/4* pi)獨(dú)立變量可用如下形式指定： dsolve( Dy=1+y2 soluti on tody/dvans=ta n(v+1/4* pi)讓我們舉一個(gè)二階微分方程的例子，該方程有兩個(gè)初始條件:2=cos(2x)-y 業(yè)(0)=0 y(0)=1 dx2dx y=dsolve( D

16、2y=cos(2*x)-y ， Dy(0)=0 ， y(0)=1 ) y=-2/3*cos(x)A2+1/3+4/3*cos(x) y=si mp le(y) % y looks like it can be simp lified y=-1/3*cos(2*x)+4/3*cos(x)通常，要求解的微分方程含有一階以上的項(xiàng)，并以下述的形式表示:烏-2 業(yè)-3y=0dx dx通解為: y=dsolve( D2y-2Dy-3*y=0 ) y=C1*ex p(-x)+C2*ex p(3*x)加上初始條件：y(0)=0和y(i)=i可得到：,y(o)=o , y(i)=i ) y=solve( D2y

17、-2Dy-3*y=0 y=1/(ex p( -1)-ex p(3)*ex p(-x)-1/(ex p( -1)-ex p( 3)*ex p(3*x) y=si mp le(y)% this looks like a can didateforsim plificatio ny=-(ex p(-x)-ex p(3*x)/(ex p(3)-ex p(-1) P retty(y) % p retty it up exp(-x)-ex p(3 x)exp(3) -exp(-l)現(xiàn)在來繪制感興趣的區(qū)域內(nèi)的結(jié)果。 ezplot(y ，-6 2)例:假設(shè)有以下三個(gè)一階常微分方程式和其初始條件y =3x2,

18、y(2)=0.5 y =2x cos(y) 2, y(0)=0.25 CfJ y =3y+ex p( 2x), y(0)=3對應(yīng)上述常微分方程式的符號運(yùn)算式為:soln_1 = dsolve(Dy = 3*x2,y(2)=0.5) ans=XA3-7.500000000000000ez plot(soln_1,2,4) %看看這個(gè)函數(shù)的長相soln_2 = dsolve(Dy = 2*x*cos(y)2,y(0) = pi/4) ans=ata n(xA2+1) soln_3 = dsolve(Dy = 3*y + exp (2*x), y(0) = 3) ans= -ex p(2*x)+4*

19、ex p(3*x)微分方程組函數(shù)dsolve也可同時(shí)處理若干個(gè)微分方程式，下面有兩個(gè)線性一階方程。f =3f+4gdx包=-4f+3g dx通解為: f，g=dsolve ( Df=3*f+4*g ，Dg=-4*f+3*g )f=C1*ex p(3*x)*si n(4*x)+C2*ex p(3*x)*cos(4*x)g=-C2*ex p(3*x)*si n(4*x)+C1*ex p(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， g(0)=1 )f=ex p(3*x

20、)*si n(4*x)g=ex p(3*x)*cos(4*x)微分和積分微分和積分是微積分學(xué)研究和應(yīng)用的核心，并廣泛地用在許多工程學(xué)科。MATLA符號工具能幫助解決許多這類問題。微分符號表達(dá)式的微分以四種形式利用函數(shù) diff ： f= a*x3+x2-b*x-cexpressionf=a*x3+x2-b*x-c% define a symbolic diff (f) % differentiate with respect to the default variable xans=3*a*xA2+2*x-b diff(f ， a ) % differentiate with respect

21、to aans=x3 diff(f ， 2) % differentiate twice with respect to xans=6*a*x+2 diff(fto aans=0a ，2) % differentiatetwice with respect函數(shù)diff也可對數(shù)組進(jìn)行運(yùn)算。如果F是符號向量或數(shù)組，diff(F)對數(shù)組內(nèi) 的各個(gè)元素進(jìn)行微分。 F=sym( a*x symbolic arrayF= a*x c*x3 ，b*x2 ；c*x3d*s ) % create a diff(F) %,b*x2d*sdifferentiate the element with respect

22、to xans= a，2*b*x03*c*x2，注意函數(shù)diff也用在MATLAB計(jì)算數(shù)值向量或矩陣的數(shù)值差分。對于一個(gè)數(shù) 值向量或矩陣M diff(M)計(jì)算M(2: m,: )-M(1 : m-1，:)的數(shù)值差分，如下 所示： m=(1 : 8)42) % create a vectorM=1 4 9 16 25 36 49 64 diff(M) % find the differences between elements ans=3 5 7 9 11 13 15如果diff的表達(dá)式或可變參量是數(shù)值，MATLA就非常巧妙地計(jì)算其數(shù)值差分；如果參量是符號字符串或變量，MATLA對其表達(dá)式進(jìn)行

23、微分。積分積分函數(shù)int(f)，其中f是一符號表達(dá)式，它力圖求出另一符號表達(dá)式F使diff(F)=f 。正如從研究微分學(xué)所了解的，積分比微分復(fù)雜得多。積分或逆求導(dǎo) 不一定是以封閉形式存在， 或許存在但軟件也許找不到， 或者軟件可明顯地求解， 但超過內(nèi)存或時(shí)間限制。當(dāng)MATLA不能找到逆導(dǎo)數(shù)時(shí)，它將返回未經(jīng)計(jì)算的命令。 in t( log(x)/ex 9(x2) ) % atte mpt to i ntegrate ans=log(x)/ex p(xA2)同微分一樣，積分函數(shù)有多種形式。形式 int(f) 相對于缺省的獨(dú)立變量求逆 導(dǎo)數(shù)；形式(f， s )相對于符號變量S積分；形式int(f ，

24、a，b)和int(f ， s，a，b)，a，b是數(shù)值，求解符號表達(dá)式從a到b的定積分；形式int(f ， m n )和形式int(f ， s ， m ， n )，其中m, n是符號變量，求解符號表達(dá)式從m到n的定積分。 f= sin(s+2*x) % crate a symbolic function f=sin(s+2*x) int(f) % integrate with respect to x ans=-1/2*cos(s+2*x) int(f ， s ) % integrate with respect to s ans=-cos(s+2*x) int(f ，pi/2 ，pi) % i

25、ntegrate with respect to x from兀/2 to兀ans=-cos(x) int(f， s ，pi/2 ，pi) % integrate with respect to s from 兀 /2 to 兀ans=cos(2*x)-si n(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)正如函數(shù)diff 樣，積分函數(shù)int對符號數(shù)組的每一個(gè)元素進(jìn)行運(yùn)算。 F=sym( a*x， b*x2 ； c*x3，d*s ) % crea

26、te a symbolic arrayF=a*x ，b*x2c*x3 ，d*s diff(F) % ubtegrate the array eleme nts with res pect to x ans=1/2*a*x2，1/3*b*xA31/4*c*x4， d*s*xdiff函數(shù)用以演算一函數(shù)的微分項(xiàng)，相關(guān)的函數(shù)語法有下列4個(gè):diff(f)傳回f對預(yù)設(shè)獨(dú)立變數(shù)的一次微分值傳回f對獨(dú)立變數(shù)t的一次微分值diff(f, n)傳回f對預(yù)設(shè)獨(dú)立變數(shù)的n次微分值diff(f,t,n)傳回 f 對獨(dú)立變數(shù) t 的 n 次微分值先定義下列三個(gè)方程式，接著再演算其微分項(xiàng)：S1 = 6*xA3-4*xA2

27、+b*x-5;S2 = sin(a);S3 = (1 -護(hù)3”(1 + 護(hù)4);diff(S1) ans= 18*x2-8*x+b diff(S1,2) ans= 36*x-8 diff(S1,b) ans=diff(S2) ans= cos(a) diff(S3) ans= -3*護(hù)2心+護(hù)4)-4*(1-護(hù)3”(1+護(hù)4)八2*護(hù)3 simplify(diff(S3) ans=護(hù)2*(-3+護(hù)4-4*圳(1+護(hù)4)八2這個(gè)函數(shù)要找出一符號式 F 使得不存在的int 函數(shù)用以演算一函數(shù)的積分項(xiàng)， diff(F)=f 。如果積分式的解析式 (analytical form, closed fo

28、rm) 話或是MATLA無法找到，則int傳回原輸入的符號式。相關(guān)的函數(shù)語法有下列 4 個(gè)：int(f) 傳回 f 對預(yù)設(shè)獨(dú)立變數(shù)的積分值int(f,t)傳回 f 對獨(dú)立變數(shù) t 的積分值int(f,a,b)式傳回f對預(yù)設(shè)獨(dú)立變數(shù)的積分值，積分區(qū)間為a,b , a和b為數(shù)值int(f,t,a,b)值式傳回f對獨(dú)立變數(shù)t的積分值，積分區(qū)間為a,b , a和b為數(shù)int(f,m,n)式傳回f對預(yù)設(shè)變數(shù)的積分值，積分區(qū)間為m,n，m和n為符號我們示范幾個(gè)例子：S1 = 6*x3-4*x2+b*x-5;S2 = sin(a);S3 = sqrt(x);int(S1)ans=3/2*x4-4/3*x3+

29、1/2*b*x2-5*xint(S2)ans=-cos(a)int(S3)ans=2/3*xA(3/2) i nt(S3,a,b) ans= 2/3*bA(3/2)- 2/3*aA(3/2) i nt(S3,0.5,0.6) ans= 2/25*15八(1/2)-1/6*2八(1/2)numeric(int(S3,0.5,0.6)%使用numeric函數(shù)可以計(jì)算積分的數(shù)值 ans= 0.074 1數(shù)值積分先考慮一個(gè)積分式的數(shù)學(xué)式如下:其中a, b分別為這個(gè)積分式的上限及下限，f(x)則為要積分的函樹。要求解上述 的積分式，必須設(shè)定a, b和f(x)。以MATLAB的積分函數(shù)求解的過程亦同， 也要定義f(x)及設(shè)定a,b，還須設(shè)定在區(qū)間a,b之間離散 點(diǎn)(discretized points) 數(shù)目，剩下的工作就是選擇精度不同的積分法來求解了。梯形法MATLAB提供最簡單的積分函數(shù)是梯形法trapz,我們先說明梯形法語法 trapz(x,y)，其中x,y分別代