MATLAB解方程與函數(shù)極值課件

1、第4章 MATLAB解方程與函數(shù)極值4.1 線性方程組求解4.2 非線性方程數(shù)值求解4.3 常微分方程初值問題的數(shù)值解法4.4 無約束優(yōu)化問題4.5 約束優(yōu)化問題10/10/20221第4章 MATLAB解方程與函數(shù)極值10/9/202214.1 線性方程組求解4.1.1 直接解法 1利用左除運算符的直接解法 對于線性方程組Ax=b，可以利用左除運算符“”求解： x=Ab 對于線性方程組xA=b，可以利用右除運算符“/”求解： x=A/b10/10/202224.1 線性方程組求解10/9/20222例4-1：用直接解法求解下列線性方程組。命令如下：A=2,1,-5,1;1,-5,0,7;0,

2、2,1,-1;1,6,-1,-4;b=13,-9,6,0;x=Ab結(jié)果：x = -66.5556 25.6667 -18.7778 26.555610/10/20223例4-1：用直接解法求解下列線性方程組。10/9/202232利用矩陣的分解求解線性方程組 矩陣分解是指根據(jù)一定的原理用某種算法將一個矩陣分解成若干個矩陣的乘積。常見的矩陣分解有LU分解、QR分解、Cholesky分解，以及Schur分解、Hessenberg分解、奇異分解等。10/10/202242利用矩陣的分解求解線性方程組10/9/20224(1) LU分解 矩陣的LU分解就是將一個矩陣表示為一個交換下三角矩陣和一個上三角

3、矩陣的乘積形式。線性代數(shù)中已經(jīng)證明，只要方陣A是非奇異的，LU分解總是可以進行的。 MATLAB提供的lu函數(shù)用于對矩陣進行LU分解，其調(diào)用格式為： L,U=lu(A)：產(chǎn)生一個上三角陣U和一個變換形式的下三角陣L(行交換)，使之滿足A=LU。注意，這里的矩陣A必須是方陣。 L,U,P=lu(A)：產(chǎn)生一個上三角陣U和一個下三角陣L以及一個置換矩陣P，使之滿足PA=LU。當然矩陣A同樣必須是方陣。 實現(xiàn)LU分解后，線性方程組Ax=b的解x=U(Lb)或x=U(LP*b)，這樣可以大大提高運算速度。10/10/20225(1) LU分解10/9/20225例：用LU分解求解例4-1中的線性方程組

4、。命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;L,U=lu(A);x=U(Lb)或采用LU分解的第2種格式，命令如下：L,U ,P=lu(A);x=U(LP*b)10/10/20226例：用LU分解求解例4-1中的線性方程組。10/9/2022 (2) QR分解 對矩陣A進行QR分解，就是把A分解為一個正交矩陣Q和一個上三角矩陣R的乘積形式。QR分解只能對方陣進行。MATLAB的函數(shù)qr可用于對矩陣進行QR分解，其調(diào)用格式為： Q,R=qr(A)：產(chǎn)生一個一個正交矩陣Q和一個上三角矩陣R，使之滿足X=QR。 Q,R,E=qr(

5、A)：產(chǎn)生一個一個正交矩陣Q、一個上三角矩陣R以及一個置換矩陣E，使之滿足AE=QR。實現(xiàn)QR分解后，線性方程組Ax=b的解x=R(Qb)或x=E(R(Qb)。10/10/20227 (2) QR分解10/9/20227例： 用QR分解求解例4-1中的線性方程組。命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;Q,R=qr(A);x=R(Qb)或采用QR分解的第2種格式，命令如下：Q,R,E=qr(A);x=E*(R(Qb)10/10/20228例： 用QR分解求解例4-1中的線性方程組。10/9/202 (3) Cholesky

6、分解 如果矩陣A是對稱正定的，則Cholesky分解將矩陣A分解成一個下三角矩陣和上三角矩陣的乘積。設(shè)上三角矩陣為R，則下三角矩陣為其轉(zhuǎn)置，即A=RR。MATLAB函數(shù)chol(A)用于對矩陣A進行Cholesky分解，其調(diào)用格式為： R=chol(A)：產(chǎn)生一個上三角陣R，使RR=A。若A為非對稱正定，則輸出一個出錯信息。 實現(xiàn)Cholesky分解后，線性方程組Ax=b變成RRx=b，所以x=R(Rb)。10/10/20229 (3) Cholesky分解10/9/20229例： 用Cholesky分解求解例4-1中的線性方程組。命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,

7、-1;1,6,-1,-4;b=13,-9,6,0;R=chol(A)? Error using = cholMatrix must be positive definite命令執(zhí)行時，出現(xiàn)錯誤信息，說明A為非正定矩陣。10/10/202210例： 用Cholesky分解求解例4-1中的線性方程組。104.1.2 迭代解法 迭代解法非常適合求解大型系數(shù)矩陣的方程組。在數(shù)值分析中，迭代解法主要包括 Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和兩步迭代法。1Jacobi迭代法 對于線性方程組Ax=b，如果A中aii0(i=1,2,n)，則可將A分解為A=D-L-U，其中D為對角

8、陣，其元素為A的對角元素，L與U為A的下三角陣和上三角陣，于是Ax=b化為： x=D-1(L+U)x+D-1b 與之對應(yīng)的迭代公式為： x(k+1)=D-1(L+U)x(k)+D-1b 這就是Jacobi迭代公式。如果序列x(k+1)收斂于x，則x必是方程Ax=b的解。 Jacobi迭代法收斂的充分必要條件是D-1(L+U)最大特征值的絕對值小于1。10/10/2022114.1.2 迭代解法10/9/202211Jacobi迭代法的MATLAB函數(shù)原文件Jacobi.m如下：function y,n=jacobi(A,b,x0,eps)if nargin=3 eps=1.0e-6;elsei

9、f nargin=eps x0=y; y=B*x0+f; n=n+1;end10/10/202212Jacobi迭代法的MATLAB函數(shù)原文件Jacobi.m如例：用Jacobi迭代法求解下列線性方程組。設(shè)迭代初值為0，迭代精度為10-6。 在命令中調(diào)用函數(shù)文件Jacobi.m，命令如下：A=10,-1,0;-1,10,-2;0,-2,10;b=9,7,6;x,n=jacobi(A,b,0,0,0,1.0e-6)結(jié)果：x = 0.9958 0.9579 0.7916n = 1110/10/202213例：用Jacobi迭代法求解下列線性方程組。設(shè)迭代初值為0，2Gauss-Serdel迭代法

10、在Jacobi迭代過程中，計算時，已經(jīng)得到，不必再用，即原來的迭代公式Dx(k+1)=(L+U)x(k)+b可以改進為Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到： x(k+1)=(D-L)-1Ux(k)+(D-L)-1b 該式即為Gauss-Serdel迭代公式。和Jacobi迭代相比，Gauss-Serdel迭代用新分量代替舊分量，精度會高些。 Gauss-Serdel迭代法收斂的充分必要條件是(D-L)-1U最大特征值的絕對值小于1。10/10/2022142Gauss-Serdel迭代法10/9/202214Gauss-Serdel迭代法的MATLAB函數(shù)原文件gausei

11、del.m如下：function y,n=gauseidel(A,b,x0,eps)if nargin=3 eps=1.0e-6;elseif nargin=eps x0=y; y=G*x0+f; n=n+1;end10/10/202215Gauss-Serdel迭代法的MATLAB函數(shù)原文件gau例：用Gauss-Serdel迭代法求解下列線性方程組。設(shè)迭代初值為0，迭代精度為10-6。 在命令中調(diào)用函數(shù)文件gauseidel.m，命令如下： A=10,-1,0;-1,10,-2;0,-2,10; b=9,7,6; x,n=gauseidel(A,b,0,0,0,1.0e-6)結(jié)果：x =

12、0.9958 0.9579 0.7916n = 710/10/202216例：用Gauss-Serdel迭代法求解下列線性方程組。設(shè)迭 若J法與GS法均收斂，則GS法比J法約快一倍，但也可能J法收斂而GS法不收斂或相反。 10/10/202217 若J法與GS法均收斂，則GS法比J法約快一例: 分別用Jacobi迭代和Gauss-Serdel迭代法求解下列線性方程組，看是否收斂。命令如下：a=1,2,-2;1,1,1;2,2,1;b=9;7;6;x,n=jacobi(a,b,0;0;0)x,n=gauseidel(a,b,0;0;0)x = -27 26 8n = 4x = NaN NaN N

13、aNn = 101210/10/202218例: 分別用Jacobi迭代和Gauss-Serdel迭代法4.2 非線性方程數(shù)值求解4.2.1 單變量非線性方程求解 在MATLAB中提供了一個fzero函數(shù)，可以用來求單變量非線性方程的根。該函數(shù)的調(diào)用格式為： z=fzero(fname,x0,tol,trace) 其中fname是待求根的函數(shù)文件名，x0為搜索的起點。一個函數(shù)可能有多個根，但fzero函數(shù)只給出離x0最近的那個根。tol控制結(jié)果的相對精度，缺省時取tol=eps，trace指定迭代信息是否在運算中顯示，為1時顯示，為0時不顯示，缺省時取trace=0。10/10/2022194

14、.2 非線性方程數(shù)值求解10/9/202219 例：求f(x)=x-10 x+2=0在x0=0.5附近的根。 步驟如下： (1) 建立函數(shù)文件funx.m。 function fx=funx(x) fx=x-10.x+2; (2) 調(diào)用fzero函數(shù)求根。 z=fzero(funx,0.5) z = 0.375810/10/202220 例：求f(x)=x-10 x+2=0在x0=0.5附近4.2.2 非線性方程組的求解 對于非線性方程組F(X)=0，用fsolve函數(shù)求其數(shù)值解。fsolve函數(shù)的調(diào)用格式為： X=fsolve(fun,X0,option) 其中X為返回的解，fun是用于定義

15、需求解的非線性方程組的函數(shù)文件名，X0是求根過程的初值，option為最優(yōu)化工具箱的選項設(shè)定。最優(yōu)化工具箱提供了20多個選項，用戶可以使用optimset命令將它們顯示出來。如果想改變其中某個選項，則可以調(diào)用optimset()函數(shù)來完成。例如，Display選項決定函數(shù)調(diào)用時中間結(jié)果的顯示方式，其中off為不顯示，iter表示每步都顯示，final只顯示最終結(jié)果。optimset(Display,off)將設(shè)定Display選項為off。10/10/2022214.2.2 非線性方程組的求解10/9/202221 例：求下列非線性方程組在(0.5,0.5) 附近的數(shù)值解。 x-0.6*sin

16、(x)-0.3*cos(y)=0; y-0.6*cos(x)+0.3*sin(y)=0; (1) 建立函數(shù)文件myfun.m。function q=myfun(x)q(1)=x(1)-0.6*sin(x(1)-0.3*cos(x(2);q(2)=x(2)-0.6*cos(x(1)+0.3*sin(x(2); (2) 在給定的初值x0=0.5,y0=0.5下，調(diào)用fsolve函數(shù)求方程的根。x=fsolve(myfun,0.5,0.5,optimset(Display,off)x = 0.6354 0.373410/10/202222 例：求下列非線性方程組在(0.5,0.5) 附近的數(shù)將求得的

17、解代回原方程，可以檢驗結(jié)果是否正確，命令如下：q=myfun(x)q = 1.0e-009 * 0.2375 0.2957 可見得到了較高精度的結(jié)果。10/10/202223將求得的解代回原方程，可以檢驗結(jié)果是否正確，命令如下：10/4.3 常微分方程初值問題的數(shù)值解法4.3.1 龍格庫塔法簡介4.3.2 龍格庫塔法的實現(xiàn) 基于龍格庫塔法，MATLAB提供了求常微分方程數(shù)值解的函數(shù)，一般調(diào)用格式為： t,y=ode23(fname,tspan,y0) t,y=ode45(fname,tspan,y0) 其中fname是定義f(t,y)的函數(shù)文件名，該函數(shù)文件必須返回一個列向量。tspan形式為

18、t0,tf,表示求解區(qū)間。y0是初始狀態(tài)列向量。t和y分別給出時間向量和相應(yīng)的狀態(tài)向量。10/10/2022244.3 常微分方程初值問題的數(shù)值解法10/9/202224例：設(shè)有初值問題y =(y2-t-2)/(4*(t+1)，y(0)=2，試求其數(shù)值解，并與精確解相比較(精確解為y(t)=sqrt(t+1)+1 )。 (1) 建立函數(shù)文件funt.m。 function dy=funt(t,y) dy=(y2-t-2)/4/(t+1); (2) 求解微分方程。 t0=0;tf=10; y0=2; t,y=ode23(funt,t0,tf,y0); %求數(shù)值解 y1=sqrt(t+1)+1;

19、%求精確解 結(jié)果為：10/10/202225例：設(shè)有初值問題y =(y2-t-2)/(4*(t+1 t = 0 0.3200 0.8200 1.3200 1.8200 2.3200 2.8200 3.3200 3.8200 4.3200 4.8200 5.0000 y =2.0000 2.1490 2.3495 2.5239 2.6803 2.8234 2.9561 3.0805 3.1978 3.3093 3.4157 3.4529 y1 =2.0000 2.1489 2.3491 2.5232 2.6793 2.8221 2.9545 3.0785 3.1954 3.3065 3.4125

20、 3.4495 y為數(shù)值解，y1為精確值，顯然兩者近似。10/10/20222610/9/202226例: 求解y-7(1-y2)y+y=0，y(0)=0.8，y (0)=0，并畫出解的圖形。 令x1=y， x2=y，x1=x2，x2=7(1-x12)x2-x1， x1(0)=0.8，x2(0)=0 建立函數(shù)文件ff.m function dx=ff(t,x) dx=x(2); 7*(1-x(1)2)*x(2)-x(1) 再編寫m文件求解 x0=0.8;0; t,x=ode45( ff ,0,40,x0); y=x(:,1); dy=x(:,2); plot(t,y,t,dy)10/10/20

21、2227例: 求解y-7(1-y2)y+y=0，y(0)=0建立函數(shù)文件ff.m function dx=ff(t,x) dx=x(2); x(3);-exp(t)-2*x(3)-3*x(2); 再編寫m文件求解 x0=0.8;0;1; t,x=ode45( ff ,0,40,x0); y=x(:,1); dy=x(:,2); d2y=x(:,3); plot(t,y,t,dy,t,d2y)10/10/202228建立函數(shù)文件ff.m10/9/202228 4.4 無約束優(yōu)化問題 MATLAB提供了基于單純形算法求解函數(shù)極值的函數(shù)fminbnd和fminunc，它們分別用于單變量函數(shù)和多變量函

22、數(shù)的無約束優(yōu)化（最小值）問題，其調(diào)用格式為： x,fval=fminbnd(fname,x1,x2) x,fval= =fminunc(fname,x0) 這兩個函數(shù)的調(diào)用格式相似。其中fminbnd函數(shù)用于求單變量函數(shù)的最小值點。fname是被最小化的目標函數(shù)名，x1和x2限定自變量的取值范圍。fminunc函數(shù)用于求多變量函數(shù)的最小值點，x0是求解的初始值向量。x返回最小值點，fval返回最小值。10/10/202229 4.4 無約束優(yōu)化問題10/9/202229 MATLAB沒有專門提供求函數(shù)最大值的函數(shù)，但只要注意到-f(x)在區(qū)間(a,b)上的最小值就是f(x)在(a,b)的最大值，所以fminbnd(-f,x1,x2)返回函數(shù)f(x)在區(qū)間(x1,x2)上的最大值。10/10/202230 MATLAB沒有專門提供求函數(shù)最 例: 求f(x)=x3-2x-5在0,5內(nèi)的最小值點以及最小值。 (1) 建立函數(shù)文件mymin.m。 function fx=mymin(x) fx=x.3-2*x-5; (2) 調(diào)用fmin