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

1、 第第7講講 MATLAB解方程解方程與函數(shù)極值與函數(shù)極值7.1 線性方程組求解線性方程組求解7.2 非線性方程數(shù)值求解非線性方程數(shù)值求解7.3 常微分方程初值問題的數(shù)值常微分方程初值問題的數(shù)值解法解法7.4 無約束最優(yōu)化問題求解無約束最優(yōu)化問題求解7.1 線性方程組求解線性方程組求解7.1.1 直接解法直接解法1利用左除運(yùn)算符的直接解法利用左除運(yùn)算符的直接解法 對于線性方程組對于線性方程組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ù)一定的原理用某種算法將一個矩陣分解成若干個矩陣的乘積。將一個矩陣分解成若干個矩陣的乘積。 常見的矩陣分解有常見的矩陣分解有LU分解分解、QR分解分解、Cholesky分解分解，以及，以及Schur分解分解、Hessenberg分解分解、奇異分解奇異分解等。下面只簡等。下面只簡單介紹前兩個和單介紹前兩個和(1) LU分解分解矩陣的矩陣的L

3、U分解就是將一個矩陣表示為一個分解就是將一個矩陣表示為一個下下三角矩陣三角矩陣L和一個上三角矩陣和一個上三角矩陣U的乘積的乘積形式。形式。A是非奇異的。其調(diào)用格式為：是非奇異的。其調(diào)用格式為：L,U=lu(A) L,U,P=lu(A)P是置換矩陣，使之滿足是置換矩陣，使之滿足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分解為分解為一個正交矩陣一個正交矩陣Q和一個上三角矩陣和一個上三角矩陣R的乘積的乘積。A是方陣進(jìn)行。其調(diào)用格式為：是方陣進(jìn)行。其調(diào)用格式為：Q,R=qr(A)Q,R,E=qr(A)E是置換矩陣，使之滿足是置換矩陣，使之滿足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為搜索的起點。為搜索的起點。tol控制結(jié)果控制結(jié)果

6、的相對精度，缺省時取的相對精度，缺省時取tol=eps (10-6)，trace 指定指定迭代信息是否在運(yùn)算中顯示，為迭代信息是否在運(yùn)算中顯示，為1時顯示，為時顯示，為0時時不顯示，缺省時取不顯示，缺省時取trace=0。注：一個函數(shù)可能有多個根，但注：一個函數(shù)可能有多個根，但fzero函數(shù)函數(shù)只給出離只給出離x0最近的那個根最近的那個根。 例例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是求根過程的初值，是求根過程的初值，option為最優(yōu)化工為最優(yōu)化工具箱的選項設(shè)定。具箱的選項設(shè)定。最優(yōu)化工具箱提供了最優(yōu)化工具箱提供了20多

8、個選項，用戶可以使用多個選項，用戶可以使用optimset命令將它們顯示出來。如果想改變其中某命令將它們顯示出來。如果想改變其中某個選項，則可以調(diào)用個選項，則可以調(diào)用optimset()函數(shù)來完成。例如，函數(shù)來完成。例如，Display選項決定函數(shù)調(diào)用時中間結(jié)果的顯示方選項決定函數(shù)調(diào)用時中間結(jié)果的顯示方式，其中式，其中off為不顯示，為不顯示，iter表示每步都顯示，表示每步都顯示，final只顯示最終結(jié)果。只顯示最終結(jié)果。optimset(Display, off)將設(shè)定將設(shè)定Display選項為選項為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、驗結(jié)果是否將求得的解代回原方程，可以檢驗結(jié)果是否正確，命令如下：正確，命令如下：q=myfun(x)q = 1.0e-009 * 0.2375 0.2957可見得到了較高精度的結(jié)果。可見得到了較高精度的結(jié)果。7.3 常微分方程初值問題的數(shù)值解法常微分方程初值問題的數(shù)值解法考慮常微分方程的初值問題考慮常微分方程的初值問題 , 所謂其數(shù)值解法，就是求所謂其數(shù)值解法，就是求它的解它的解y(t)在節(jié)點在節(jié)點t0t1tm處的近似值處的近似值y0,y1,ym的方法。的方法。所求得的所求得的y0,y1,ym稱為常微分方程初值問稱為常微分方程初值問題的數(shù)值解。一般采用等距節(jié)點題的數(shù)值解。一般采用等距節(jié)點tn=

11、t0+nh, n=0,1,m, h稱為稱為步長步長。常用的有常用的有歐拉歐拉(Euler)法法、龍格庫塔龍格庫塔(Runge- Kutta)法法、線性多步法線性多步法、預(yù)報校正法預(yù)報校正法等。等。Tttytfy0),(00)(yty7.3.1 龍格庫塔法遞推公式龍格庫塔法遞推公式其中其中)22(6)(432110kkkkhyyihtyii),()2,2()2,2(),(311421131112111hkyhtfkkhyhtfkkhyhtfkytfkiiiiiiii7.3.2 龍格庫塔法的實現(xiàn)龍格庫塔法的實現(xiàn)基于龍格庫塔法，基于龍格庫塔法，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ù)文件必須返回一個列向量。文件必須返回一個列向量。tspan形式為形式為t0,tf,表示求解區(qū)間。表示求解區(qū)間。y0是初始狀態(tài)列向量。是初始狀態(tài)列向量。t和和y分別給出時間向量和相應(yīng)的狀態(tài)向量。分別給出時間向量和相應(yīng)的狀態(tài)向量。例例7-6 設(shè)有初值問題設(shè)有初值問題 求其數(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-); %通過圖形來比較通過圖形來比較 2)0(,100 ,) 1(422ytttyy11 ty結(jié)果如下結(jié)果如下數(shù)值解為藍(lán)色圓點，精確值為紅色實線，顯數(shù)值解為藍(lán)色圓點，精確值為紅色實線，顯然兩者結(jié)果近似。然兩者結(jié)果近似。例例7-7 已知一個二階線性系統(tǒng)微分方程為已知一個

14、二階線性系統(tǒng)微分方程為 其中其中a=2，繪制系統(tǒng)的時間響應(yīng)曲線和相平面圖，繪制系統(tǒng)的時間響應(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 無約束最優(yōu)化問題求解無約束最優(yōu)化問題求解MATLAB提供了提供了3個最小值函數(shù)，調(diào)用格式為：個最小值函數(shù)，調(diào)用格式為：(1) x,fval=fminbnd(fname,x1,x2)求一元函數(shù)在區(qū)間求一元

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