第7講MATLAB解方程與函數(shù)極值

1、 第第7講講 MATLAB解方程解方程與函數(shù)極值與函數(shù)極值7.1 線性方程組求解線性方程組求解7.2 非線性方程數(shù)值求解非線性方程數(shù)值求解7.3 常微分方程初值問(wèn)題的數(shù)值常微分方程初值問(wèn)題的數(shù)值解法解法7.4 無(wú)約束最優(yōu)化問(wèn)題求解無(wú)約束最優(yōu)化問(wèn)題求解7.1 線性方程組求解線性方程組求解7.1.1 直接解法直接解法1利用左除運(yùn)算符的直接解法利用左除運(yùn)算符的直接解法 對(duì)于線性方程組對(duì)于線性方程組Ax=b， 可以利用左除運(yùn)算符可以利用左除運(yùn)算符“”求解：求解：x=Ab例例7-1 用直接解法求解下列線性方程組用直接解法求解下列線性方程組 例例命令如下：命令如下：A=2,1,-5,1;1,-5,0,7;

2、0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;x=Ab04662975135243214324214321xxxxxxxxxxxxxx2利用利用矩陣的分解矩陣的分解求解線性方程組求解線性方程組 矩陣分解是指根據(jù)一定的原理用某種算法矩陣分解是指根據(jù)一定的原理用某種算法將一個(gè)矩陣分解成若干個(gè)矩陣的乘積。將一個(gè)矩陣分解成若干個(gè)矩陣的乘積。 常見(jiàn)的矩陣分解有常見(jiàn)的矩陣分解有LU分解分解、QR分解分解、Cholesky分解分解，以及，以及Schur分解分解、Hessenberg分解分解、奇異分解奇異分解等。下面只簡(jiǎn)等。下面只簡(jiǎn)單介紹前兩個(gè)和單介紹前兩個(gè)和(1) LU分解分解矩陣的矩陣的L

3、U分解就是將一個(gè)矩陣表示為一個(gè)分解就是將一個(gè)矩陣表示為一個(gè)下下三角矩陣三角矩陣L和一個(gè)上三角矩陣和一個(gè)上三角矩陣U的乘積的乘積形式。形式。A是非奇異的。其調(diào)用格式為：是非奇異的。其調(diào)用格式為：L,U=lu(A) L,U,P=lu(A)P是置換矩陣，使之滿(mǎn)足是置換矩陣，使之滿(mǎn)足PA=LU。LU分解可以大大提高運(yùn)算速度。分解可以大大提高運(yùn)算速度。例例7-2 用用LU分解求解分解求解例例7-1中的線性方程組。中的線性方程組。命令如下：命令如下：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分解的

4、第分解的第2種格式，命令如下：種格式，命令如下：L,U ,P=lu(A);x=U(LP*b) (2) QR分解分解把把A分解為分解為一個(gè)正交矩陣一個(gè)正交矩陣Q和一個(gè)上三角矩陣和一個(gè)上三角矩陣R的乘積的乘積。A是方陣進(jìn)行。其調(diào)用格式為：是方陣進(jìn)行。其調(diào)用格式為：Q,R=qr(A)Q,R,E=qr(A)E是置換矩陣，使之滿(mǎn)足是置換矩陣，使之滿(mǎn)足AE=QR。例例7-3 用用QR分解求解例分解求解例7-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)或采用或采用

5、QR分解的第分解的第2種格式，命令如下：種格式，命令如下：Q,R,E=qr(A);x=E*(R(Qb)Ex:分別用矩陣求逆、矩陣除法以及矩陣分解分別用矩陣求逆、矩陣除法以及矩陣分解求下列線性方程組的解：求下列線性方程組的解： 57347310532zyxzyxzyx7.2 非線性方程數(shù)值求解非線性方程數(shù)值求解7.2.1 單變量非線性方程求解單變量非線性方程求解fzero函數(shù)求函數(shù)求單變量單變量非線性方程的根。調(diào)用格式為：非線性方程的根。調(diào)用格式為： z=fzero(fname,x0,tol,trace)fname是函數(shù)文件名，是函數(shù)文件名，x0為搜索的起點(diǎn)。為搜索的起點(diǎn)。tol控制結(jié)果控制結(jié)果

6、的相對(duì)精度，缺省時(shí)取的相對(duì)精度，缺省時(shí)取tol=eps (10-6)，trace 指定指定迭代信息是否在運(yùn)算中顯示，為迭代信息是否在運(yùn)算中顯示，為1時(shí)顯示，為時(shí)顯示，為0時(shí)時(shí)不顯示，缺省時(shí)取不顯示，缺省時(shí)取trace=0。注：一個(gè)函數(shù)可能有多個(gè)根，但注：一個(gè)函數(shù)可能有多個(gè)根，但fzero函數(shù)函數(shù)只給出離只給出離x0最近的那個(gè)根最近的那個(gè)根。 例例7-4 求求f(x)=x-10 x+2=0在在x0=0.5附近的根。附近的根。 步驟如下：步驟如下： (1) 建立函數(shù)文件建立函數(shù)文件funx.m。 function fx=funx(x) fx=x-10.x+2; (2) 調(diào)用調(diào)用fzero函數(shù)求根。

7、函數(shù)求根。 z=fzero(funx,0.5) z = 0.3758Ex: 求方程求方程 在在 附近的根。附近的根。 0sinxxx5 . 00 x7.2.2 非線性方程組的求解非線性方程組的求解非線性方程組非線性方程組F(X)=0用用fsolve函數(shù)求其數(shù)值解。調(diào)用函數(shù)求其數(shù)值解。調(diào)用格式為：格式為： X=fsolve(fun,X0,option)其中其中fun是用于定義需求解的非線性方程組的函數(shù)文是用于定義需求解的非線性方程組的函數(shù)文件名，件名，X0是求根過(guò)程的初值，是求根過(guò)程的初值，option為最優(yōu)化工為最優(yōu)化工具箱的選項(xiàng)設(shè)定。具箱的選項(xiàng)設(shè)定。最優(yōu)化工具箱提供了最優(yōu)化工具箱提供了20多

8、個(gè)選項(xiàng)，用戶(hù)可以使用多個(gè)選項(xiàng)，用戶(hù)可以使用optimset命令將它們顯示出來(lái)。如果想改變其中某命令將它們顯示出來(lái)。如果想改變其中某個(gè)選項(xiàng)，則可以調(diào)用個(gè)選項(xiàng)，則可以調(diào)用optimset()函數(shù)來(lái)完成。例如，函數(shù)來(lái)完成。例如，Display選項(xiàng)決定函數(shù)調(diào)用時(shí)中間結(jié)果的顯示方選項(xiàng)決定函數(shù)調(diào)用時(shí)中間結(jié)果的顯示方式，其中式，其中off為不顯示，為不顯示，iter表示每步都顯示，表示每步都顯示，final只顯示最終結(jié)果。只顯示最終結(jié)果。optimset(Display, off)將設(shè)定將設(shè)定Display選項(xiàng)為選項(xiàng)為off。例例7-5求下列非線性方程組求下列非線性方程組 在在(0.5,0.5) 附近的數(shù)值

9、解。附近的數(shù)值解。 (1) 建立函數(shù)文件建立函數(shù)文件myfun.m。function q=myfun(p)x=p(1);y=p(2);q(1)=x-0.6*sin(x)-0.3*cos(y);q(2)=y-0.6*cos(x)+0.3*sin(y); (2) 在給定的初值在給定的初值x0=0.5,y0=0.5下，調(diào)用下，調(diào)用fsolve函數(shù)求方程的根。函數(shù)求方程的根。x=fsolve(myfun,0.5,0.5,optimset(Display,off)0sin3 . 0cos6 . 00cos3 . 0sin6 . 0yxyyxxx = 0.6354 0.3734將求得的解代回原方程，可以檢

10、驗(yàn)結(jié)果是否將求得的解代回原方程，可以檢驗(yàn)結(jié)果是否正確，命令如下：正確，命令如下：q=myfun(x)q = 1.0e-009 * 0.2375 0.2957可見(jiàn)得到了較高精度的結(jié)果。可見(jiàn)得到了較高精度的結(jié)果。7.3 常微分方程初值問(wèn)題的數(shù)值解法常微分方程初值問(wèn)題的數(shù)值解法考慮常微分方程的初值問(wèn)題考慮常微分方程的初值問(wèn)題 , 所謂其數(shù)值解法，就是求所謂其數(shù)值解法，就是求它的解它的解y(t)在節(jié)點(diǎn)在節(jié)點(diǎn)t0t1tm處的近似值處的近似值y0,y1,ym的方法。的方法。所求得的所求得的y0,y1,ym稱(chēng)為常微分方程初值問(wèn)稱(chēng)為常微分方程初值問(wèn)題的數(shù)值解。一般采用等距節(jié)點(diǎn)題的數(shù)值解。一般采用等距節(jié)點(diǎn)tn=

11、t0+nh, n=0,1,m, h稱(chēng)為稱(chēng)為步長(zhǎng)步長(zhǎng)。常用的有常用的有歐拉歐拉(Euler)法法、龍格庫(kù)塔龍格庫(kù)塔(Runge- Kutta)法法、線性多步法線性多步法、預(yù)報(bào)校正法預(yù)報(bào)校正法等。等。Tttytfy0),(00)(yty7.3.1 龍格庫(kù)塔法遞推公式龍格庫(kù)塔法遞推公式其中其中)22(6)(432110kkkkhyyihtyii),()2,2()2,2(),(311421131112111hkyhtfkkhyhtfkkhyhtfkytfkiiiiiiii7.3.2 龍格庫(kù)塔法的實(shí)現(xiàn)龍格庫(kù)塔法的實(shí)現(xiàn)基于龍格庫(kù)塔法，基于龍格庫(kù)塔法，MATLAB提供了求常微提供了求常微分方程數(shù)值解的函數(shù)，

12、一般調(diào)用格式為：分方程數(shù)值解的函數(shù)，一般調(diào)用格式為：t,y=ode23(fname,tspan,y0)t,y=ode45(fname,tspan,y0)其中其中fname是定義是定義f(t,y)的函數(shù)文件名，該函數(shù)的函數(shù)文件名，該函數(shù)文件必須返回一個(gè)列向量。文件必須返回一個(gè)列向量。tspan形式為形式為t0,tf,表示求解區(qū)間。表示求解區(qū)間。y0是初始狀態(tài)列向量。是初始狀態(tài)列向量。t和和y分別給出時(shí)間向量和相應(yīng)的狀態(tài)向量。分別給出時(shí)間向量和相應(yīng)的狀態(tài)向量。例例7-6 設(shè)有初值問(wèn)題設(shè)有初值問(wèn)題 求其數(shù)值解，并與精確解比較求其數(shù)值解，并與精確解比較(精確解為精確解為 )。(1) 建立函數(shù)文件建立函

13、數(shù)文件funt.m。function yp=funt(t,y)yp=(y2-t-2)/4/(t+1);(2) 求解微分方程。求解微分方程。t0=0;tf=10;y0=2;t,y=ode23(funt,t0,tf,y0); %求數(shù)值解求數(shù)值解y1=sqrt(t+1)+1; %求精確解求精確解plot(t,y, b.,t,y1, r-); %通過(guò)圖形來(lái)比較通過(guò)圖形來(lái)比較 2)0(,100 ,) 1(422ytttyy11 ty結(jié)果如下結(jié)果如下數(shù)值解為藍(lán)色圓點(diǎn)，精確值為紅色實(shí)線，顯數(shù)值解為藍(lán)色圓點(diǎn)，精確值為紅色實(shí)線，顯然兩者結(jié)果近似。然兩者結(jié)果近似。例例7-7 已知一個(gè)二階線性系統(tǒng)微分方程為已知一個(gè)

14、二階線性系統(tǒng)微分方程為 其中其中a=2，繪制系統(tǒng)的時(shí)間響應(yīng)曲線和相平面圖，繪制系統(tǒng)的時(shí)間響應(yīng)曲線和相平面圖.分析：令分析：令x2=x， x1=x，則系統(tǒng)的狀態(tài)方程：，則系統(tǒng)的狀態(tài)方程：建立函數(shù)文件：建立函數(shù)文件：1)0(, 0)0(0, 022xxaaxdtxd1)0(, 0)0(122112xxaxxxxfunction xdot=sys(t,x)xdot=-2*x(2);x(1);取t0=0，tf=20,求微分方程的解：t0=0;tf=20;t,x=ode45(sys,t0,tf,1,0) ;t,xans = 0 1.0000 0 0.0001 1.0000 0.0001 19.9170

15、-0.9937 0.0738 20.0000 -0.9991 -0.0090subplot(1,2,1);plot(t,x(:,2); %解的曲線解的曲線t-xsubplot(1,2,2);plot(x(:,2),x(:,1); %相平面曲線相平面曲線x-xaxis equal05101520-0.8-0.6-0.4-0.200.20.40.60.8-0.500.5-1-0.500.517.4 無(wú)約束最優(yōu)化問(wèn)題求解無(wú)約束最優(yōu)化問(wèn)題求解MATLAB提供了提供了3個(gè)最小值函數(shù)，調(diào)用格式為：個(gè)最小值函數(shù)，調(diào)用格式為：(1) x,fval=fminbnd(fname,x1,x2)求一元函數(shù)在區(qū)間求一元

16、函數(shù)在區(qū)間(x1,x2)中的最小值點(diǎn)中的最小值點(diǎn)x和最小值和最小值fval.(2) x=fminsearch(fname,x1,x2)(3) x=fminunc(fname,x0)這兩個(gè)函數(shù)的調(diào)用格式相似。其中這兩個(gè)函數(shù)的調(diào)用格式相似。其中fminsearch函數(shù)函數(shù)是基于單純形算法，是基于單純形算法，fminunc是基于擬牛頓法求是基于擬牛頓法求多元函數(shù)的最小值點(diǎn)，多元函數(shù)的最小值點(diǎn)，x0是求解的初始值向量。是求解的初始值向量。注意：注意：1. 當(dāng)目標(biāo)函數(shù)的階數(shù)大于當(dāng)目標(biāo)函數(shù)的階數(shù)大于2時(shí)，使用時(shí)，使用fminunc比比f(wàn)minsearch更有效。但當(dāng)目標(biāo)函數(shù)高度更有效。但當(dāng)目標(biāo)函數(shù)高度不連續(xù)時(shí)，使用不連續(xù)時(shí)，使用fminsearch效果更好效果更好2. MATLAB沒(méi)有專(zhuān)門(mén)提供求函數(shù)最大值的函沒(méi)有專(zhuān)門(mén)提供求函數(shù)最大值的函數(shù)，但只要注意到數(shù)，但只要注意到-f(x)在區(qū)間在區(qū)間(a,b)上的最上的最小值就是小值就是f(x)在在(a,b)的最大值，所以的最大值，所以fminunc(-f,x1,x2) 返回函數(shù)返回函數(shù)f(x)在區(qū)間在區(qū)間(x1,