第8章 MATLAB方程求解

第8章MATLAB方程求解

8.1線性方程組求解8.2非線性方程數(shù)值求解8.3最優(yōu)化問題求解8.4常微分方程初值問題的數(shù)值求解第8章MATLAB方程求解8.1線性方程組求解在MATLAB中，關(guān)于線性方程組的解法一般分為兩類：一類是直接法，就是在沒有舍入誤差的情況下，通過有限步的矩陣初等運算來求得方程組的解；另一類是迭代法，就是先給定一個解的初始值，然后按照一定的迭代算法進行逐步逼近，求出更精確的近似解。第8章MATLAB方程求解8.1.1線性方程組的直接解法在MATLAB中，這些算法已經(jīng)被編成了現(xiàn)成的庫函數(shù)或運算符，因此，只需調(diào)用相應(yīng)的函數(shù)或運算符即可完成線性方程組的求解。1．利用左除運算符的直接解法線性方程組求解最簡單的方法就是使用左除運算符“\”，系統(tǒng)會自動根據(jù)輸入的系數(shù)矩陣判斷選用哪種方法進行求解。對于線性方程組Ax=b，可以利用左除運算符“\”求解：x=A\b。第8章MATLAB方程求解例8-1用直接解法求解下列線性方程組。A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';x=A\b程序運行結(jié)果為：x=-66.555625.6667-18.7778

26.5556第8章MATLAB方程求解2．利用矩陣的分解求解線性方程組（1）LU分解MATLAB提供的lu函數(shù)用于對矩陣進行LU分解，其調(diào)用格式為：①[L,U]=lu(X)：產(chǎn)生一個上三角陣U和一個變換形式的下三角陣L（行交換），使之滿足X=LU。注意，這里的矩陣X必須是方陣。②[L,U,P]=lu(X)：產(chǎn)生一個上三角陣U和一個下三角陣L以及一個置換矩陣P，使之滿足PX=LU。當(dāng)然矩陣X同樣必須是方陣。當(dāng)使用第一種格式時，矩陣L往往不是一個下三角陣，但可以通過行交換成為一個下三角陣。第8章MATLAB方程求解設(shè)：則對矩陣A進行LU分解的命令如下：>>A=[1,-1,1;5,-4,3;2,1,1];>>[L,U]=lu(A)L=0.2000-0.07691.00001.0000000.40001.00000U=5.0000-4.00003.000002.6000-0.2000000.3846為檢驗結(jié)果是否正確，輸入命令：>>LU=L*ULU=

1-115-43211第8章MATLAB方程求解利用第二種格式對矩陣A進行LU分解。>>[L,U,P]=lu(A)L=1.0000000.40001.000000.2000-0.07691.0000U=5.0000-4.00003.000002.6000-0.2000000.3846P=010001100>>LU=L*U %這種分解其乘積不為ALU=5-432111-11>>inv(P)*L*U %考慮矩陣P后的結(jié)果ans=1-115-43211實現(xiàn)LU分解后，線性方程組Ax=b的解x=U\(L\b)或x=U\(L\Pb)，這樣可以大大提高運算速度。第8章MATLAB方程求解例8-2用LU分解求解例8-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\(L\b)程序運行結(jié)果為：x=-66.555625.6667-18.777826.5556或采用LU分解的第二種格式，命令如下：>>[L,U,P]=lu(A);>>x=U\(L\P*b)第8章MATLAB方程求解（2）QR分解對矩陣X進行QR分解，就是把X分解為一個正交矩陣Q和一個上三角陣R的乘積形式。QR分解只能對方陣進行。MATLAB的函數(shù)qr可用于對矩陣進行QR分解，其調(diào)用格式為：①[Q,R]=qr(X)：產(chǎn)生一個正交矩陣Q和一個上三角陣R，使之滿足X=QR。②[Q,R,E]=qr(X)：產(chǎn)生一個正交矩陣Q、一個上三角陣R以及一個置換矩陣E，使之滿足XE=QR。第8章MATLAB方程求解設(shè)：則對矩陣A進行QR分解的命令如下：>>A=[1,-1,1;5,-4,3;2,7,10];>>[Q,R]=qr(A)Q=-0.1826-0.0956-0.9785-0.9129-0.35320.2048-0.36510.9307-0.0228R=-5.47721.2780-6.572708.02298.151700-0.5917>>QR=Q*RQR=

1.0000-1.00001.00005.0000-4.00003.00002.00007.000010.0000第8章MATLAB方程求解利用第二種格式對矩陣A進行QR分解：>>[Q,R,E]=qr(A)Q=-0.0953-0.2514-0.9632-0.2860-0.91990.2684-0.95350.30110.0158R=-10.4881-5.4347-3.432506.0385-4.2485000.4105E=001010100>>Q*R/E%驗證A=Q*R*inv(E)ans=1.0000-1.00001.00005.0000-4.00003.00002.00007.000010.0000實現(xiàn)QR分解后，線性方程組Ax=b的解x=R\(Q\b)或x=E(R\(Q\b))。第8章MATLAB方程求解例8-3用QR分解求解例8-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\(Q\b)程序運行結(jié)果為：x=-66.555625.6667-18.777826.5556或采用QR分解的第二種格式，命令如下：>>[Q,R,E]=qr(A);>>x=E*(R\(Q\b))將得到與上面同樣的結(jié)果。第8章MATLAB方程求解（3）Cholesky分解如果矩陣X是對稱正定的，則Cholesky分解將矩陣X分解成一個下三角陣和一個上三角陣的乘積。設(shè)上三角陣為R，則下三角陣為其轉(zhuǎn)置，即X=R'R。MATLAB函數(shù)chol(X)用于對矩陣X進行Cholesky分解，其調(diào)用格式為：①R=chol(X)：產(chǎn)生一個上三角陣R，使R'R=X。若X為非對稱正定，則輸出一個出錯信息。②[R,p]=chol(X)：這個命令格式將不輸出出錯信息。若X為對稱正定的，則p=0，R與上述格式得到的結(jié)果相同，否則p為一個正整數(shù)。如果X為滿秩矩陣，則R為一個階數(shù)為q=p-1的上三角陣，且滿足R‘R=X(1:q,1:q)。第8章MATLAB方程求解設(shè)：則對矩陣A進行Cholesky分解的命令如下：>>A=[2,1,1;1,2,-1;1,-1,3];>>R=chol(A)R=1.41420.70710.707101.2247-1.2247001.0000可以驗證R'R=A，命令如下：>>R'*Rans=2.00001.00001.00001.00002.0000-1.00001.0000-1.00003.0000第8章MATLAB方程求解利用第二種格式對矩陣A進行Cholesky分解：>>[R,p]=chol(A)R=1.41420.70710.707101.2247-1.2247001.0000p=0結(jié)果中p=0，這表示矩陣A是一個正定矩陣。如果試圖對一個非正定矩陣進行Cholesky分解，則將得出錯誤信息，所以，chol函數(shù)還可以用來判定矩陣是否為正定矩陣。實現(xiàn)Cholesky分解后，線性方程組Ax=b變成R'Rx=b，所以x=R\(R'\b)。第8章MATLAB方程求解例8-4用Cholesky分解求解例8-1中的線性方程組。命令如下：>>A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];>>b=[13,-9,6,0]';>>R=chol(A)錯誤使用chol矩陣必須為正定矩陣。命令執(zhí)行時，出現(xiàn)錯誤信息，說明A為非正定矩陣。第8章MATLAB方程求解8.1.2線性方程組的迭代解法迭代法是一種不斷用變量的舊值遞推新值的過程，是用計算機解決問題的一種基本方法。它利用計算機運算速度快、適合做重復(fù)性操作的特點，讓一組指令重復(fù)執(zhí)行，在每次執(zhí)行這組指令時，都從變量的原值推出它的新值。迭代解法非常適合求解大型稀疏矩陣的方程組。在數(shù)值分析中，迭代解法主要包括Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和兩步迭代法。第8章MATLAB方程求解迭代法的思想:將方程改寫為：這種形式的好處是將一組x代入右端，可以立即得到另一組x。如果兩組x相等，那么它就是方程組的解，不等時可以繼續(xù)迭代?？梢詷?gòu)造方程的迭代公式為：第8章MATLAB方程求解1．Jacobi迭代法對于線性方程組Ax=b，如果A為非奇異方陣，且aii≠0（i=1，2，…，n），則可將A分解為A=D-L-U，其中D為對角陣，其元素為A的對角元素，L與U為A的下三角陣取反和上三角陣取反：于是Ax=b轉(zhuǎn)化為：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的解。第8章MATLAB方程求解Jacobi迭代法的MATLAB函數(shù)文件jacobi.m如下：function[y,n]=jacobi(A,b,x0,ep)ifnargin==3ep=1.0e-6;elseifnargin<3errorreturnendD=diag(diag(A)); %求A的對角矩陣L=-tril(A,-1); %求A的下三角陣U=-triu(A,1); %求A的上三角陣B=D\(L+U);f=D\b;y=B*x0+f;n=1; %迭代次數(shù)whilenorm(y-x0)>=epx0=y;y=B*x0+f;n=n+1;end第8章MATLAB方程求解例8-5用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.99580.95790.7916n=11第8章MATLAB方程求解2．Gauss-Serdel迭代法在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迭代用新分量代替舊分量，精度會高些。第8章MATLAB方程求解Gauss-Serdel迭代法的MATLAB函數(shù)文件gauseidel.m如下：function[y,n]=gauseidel(A,b,x0,ep)ifnargin==3ep=1.0e-6;elseifnargin<3errorreturnendD=diag(diag(A)); %求A的對角矩陣L=-tril(A,-1); %求A的下三角陣U=-triu(A,1); %求A的上三角陣G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;%迭代次數(shù)whilenorm(y-x0)>=epx0=y;y=G*x0+f;n=n+1;end第8章MATLAB方程求解例8-6用Gauss-Serdel迭代法求解例8-5中的線性方程組。在程序中調(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=0.99580.95790.7916n=7由此可見，一般情況下Gauss-Serdel迭代比Jacobi迭代要收斂快一些。但這也不是絕對的，在某些情況下，Jacobi迭代收斂而Gauss-Serdel迭代卻可能不收斂。

第8章MATLAB方程求解例8-7分別用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=-27268n=4>>[x,n]=gauseidel(a,b,[0;0;0])x=NaNNaNNaNn=

1012可見對此方程，用Jacobi迭代收斂，而Gauss-Serdel迭代不收斂。因此，在使用迭代法時，要考慮算法的收斂性。第8章MATLAB方程求解8.1.3求線性方程組的通解①當(dāng)系數(shù)矩陣A是一個滿秩方陣時，方程Ax=b稱為恰定方程，方程有唯一解x=A-1b，這是最基本的一種情況。一般用x=A\b求解速度更快。②當(dāng)方程組右端向量b=0時，方程稱為齊次方程組。齊次方程組總有零解，因此稱解x=0為平凡解。當(dāng)系數(shù)矩陣A的秩小于n（n為方程組中未知變量的個數(shù)）時，齊次方程組有無窮多個非平凡解，其通解中包含n-rank(A)個線性無關(guān)的解向量，用MATLAB的函數(shù)null(A,'r')可求得基礎(chǔ)解系。③當(dāng)方程組右端向量b≠0時，系數(shù)矩陣的秩rank(A)與其增廣矩陣的秩rank([A,b])是判斷其是否有解的基本條件：（a）當(dāng)rank(A)=rank([A,b])=n時，方程組有唯一解：x=A\b或x=inv(A)*b。（b）當(dāng)rank(A)=rank([A,b])<n時，方程組有無窮多個解，其通解=方程組的一個特解+對應(yīng)的齊次方程組Ax=0的通解?？梢杂肁\b求得方程組的一個特解，用null(A,'r')求得該方程組所對應(yīng)的齊次方程組的基礎(chǔ)解系，基礎(chǔ)解系中包含n-rank(A)個線性無關(guān)的解向量。（c）當(dāng)rank(A)<rank([A,b])時，方程組無解。第8章MATLAB方程求解下面設(shè)計一個求解線性方程組的函數(shù)文件line_solution.m。function[x,y]=line_solution(A,b)[m,n]=size(A);y=[];ifnorm(b)>0%非齊次方程組

ifrank(A)==rank([A,b])ifrank(A)==n%有唯一解

disp('原方程組有唯一解x')x=A\b;else%方程組有無窮多個解，基礎(chǔ)解系

disp('原方程組有無窮個解，特解為x，其齊次方程組的基礎(chǔ)解系為y')x=A\b;y=null(A,'r');

endelse

disp('方程組無解')%方程組無解

x=[];endelse%齊次方程組

disp('原方程組有零解x')x=zeros(n,1);%0解

ifrank(A)<ndisp('方程組有無窮個解，基礎(chǔ)解系為y')%非0解

y=null(A,'r');endend第8章MATLAB方程求解例8-8求解方程組。程序如下：A=[1,-2,3,-1;3,-1,5,-3;2,1,2,-2];b=[1;2;3];[x,y]=line_solution(A,b)程序運行結(jié)果為：方程組無解x=[]y=[]說明該方程組無解。第8章MATLAB方程求解例8-9求方程組的通解。程序如下：formatrat%指定有理式格式輸出A=[1,1,-3,-1;3,-1,-3,4;1,5,-9,-8];b=[1,4,0]';[x,y]=line_solution(A,b);x,yformatshort%恢復(fù)默認的短格式輸出程序運行結(jié)果為：原方程組有無窮個解，特解為x，其齊次方程組的基礎(chǔ)解系為y>Inline_solution(line11)Inaaa(line4)警告:秩不足，秩=2，tol=8.837264e-15。第8章MATLAB方程求解x=00-8/153/5y=3/2-3/43/27/41001所以原方程組的通解為：X=k1

+k2+，其中k1、k2為任意常數(shù)。第8章MATLAB方程求解8.2非線性方程數(shù)值求解非線性方程的求根方法很多，常用的有牛頓迭代法，但該方法需要求原方程的導(dǎo)數(shù)，而在實際運算中這一條件有時是不能滿足的，所以又出現(xiàn)了弦截法、二分法等其他方法。在MATLAB中，非線性方程的求解和最優(yōu)化問題往往需要調(diào)用最優(yōu)化工具箱來解決。優(yōu)化工具箱提供了一系列的優(yōu)化算法函數(shù)，可用于解決工程中的最優(yōu)化問題，包括非線性方程求解、極小值問題、最小二乘問題等。第8章MATLAB方程求解8.2.1單變量非線性方程求解在MATLAB中提供了一個fzero函數(shù)，可以用來求單變量非線性方程的根。該函數(shù)的調(diào)用格式為：z=fzero(filename,x0,tol,trace)其中，filename是待求根的函數(shù)文件名，x0為搜索的起點。一個函數(shù)可能有多個根，但fzero函數(shù)只給出離x0最近的那個根。tol控制結(jié)果的相對精度，默認時取tol=eps，trace指定迭代信息是否在運算中顯示，為1時顯示，為0時不顯示，默認時取trace=0。第8章MATLAB方程求解例8-10求在x0=-5和x0=1作為迭代初值時的根。先建立函數(shù)文件fz.m。functionf=fz(x)f=x-1./x+5;然后調(diào)用fzero函數(shù)求根，命令如下：>>fzero(@fz,-5)%以-5作為迭代初值ans=-5.1926>>fzero(@fz,1)%以1作為迭代初值ans=0.1926可以繪制f(x)的圖像（如下圖所示），從中可以看出f(x)在x=-5和x=1附近的零點。第8章MATLAB方程求解第8章MATLAB方程求解8.2.2非線性方程組的求解在MATLAB的最優(yōu)化工具箱中提供了非線性方程組的求解函數(shù)fsolve，其調(diào)用格式為：X=fsolve(filename,X0,option)其中，X為返回的解，filename是用于定義需求解的非線性方程組的函數(shù)文件名，X0是求解過程的初值，option用于設(shè)置優(yōu)化工具箱的優(yōu)化參數(shù)。優(yōu)化工具箱提供了許多優(yōu)化參數(shù)選項，用戶在命令行窗口輸入下列命令，可以將優(yōu)化參數(shù)全部顯示出來：>>optimset如果希望得到某個優(yōu)化函數(shù)（例如fsolve函數(shù)）當(dāng)前的默認參數(shù)值，則可在命令行窗口輸入命令：>>optimsetfsolve如果想改變其中某個參數(shù)選項，則可以調(diào)用optimset函數(shù)來完成。例如，Display參數(shù)選項決定函數(shù)調(diào)用時中間結(jié)果的顯示方式，其中'off'為不顯示，'iter'表示每步都顯示，'final'只顯示最終結(jié)果。如果要將設(shè)定Display選項為off，可以使用命令：>>option=optimset('Display','off')。第8章MATLAB方程求解例8-11求下列方程組在（1，1，1）附近的解并對結(jié)果進行驗證。首先建立函數(shù)文件myfun.m。functionF=myfun(X)x=X(1);y=X(2);z=X(3);F(1)=sin(x)+y+z^2*exp(x);F(2)=x+y+z;F(3)=x*y*z;在給定的初值x0=1，y0=1，z0=1下，調(diào)用fsolve函數(shù)求方程的根，命令如下：>>option=optimset('Display','off');>>X=fsolve(@myfun,[1,1,1],option)X=0.0224-0.0224-0.0000將求得的解代回原方程，可以檢驗結(jié)果是否正確，命令如下：>>q=myfun(X)q=1.0e-06*-0.5931-0.00000.0006第8章MATLAB方程求解例8-12求圓和直線的兩個交點。圓：直線：該問題即為求解方程組：使用fsolve函數(shù)求解方程組時，必須先估計出方程組的根的大致范圍。所給直線的方向數(shù)是(-12，24，-14)，故其與球心在坐標原點的球面的交點大致是(-1，1，-1)和(1，-1，1)。以這兩點作為迭代初值。第8章MATLAB方程求解先建立方程組函數(shù)文件fxyz.m。functionF=fxyz(X)x=X(1);y=X(2);z=X(3);F(1)=x^2+y^2+z^2-9;F(2)=3*x+5*y+6*z;F(3)=x-3*y-6*z-1;再在MATLAB命令行窗口，輸入如下命令：>>X1=fsolve('fxyz',[-1,1,-1],optimset('Display','off'))%求第一個交點X1=-0.95082.4016-1.5259>>X2=fsolve('fxyz',[1,-1,1],optimset('Display','off'))%求第二個交點X2=1.4180-2.33611.2377第8章MATLAB方程求解8.3最優(yōu)化問題求解最優(yōu)化方法包含有多個分支，如線性規(guī)劃、整數(shù)規(guī)劃、非線性規(guī)劃、動態(tài)規(guī)劃、多目標規(guī)劃等。利用MATLAB的優(yōu)化工具箱，可以求解線性規(guī)劃、非線性規(guī)劃和多目標規(guī)劃問題。MATLAB還提供了求非線性函數(shù)最小值問題的求解方法，為優(yōu)化方法在工程中的實際應(yīng)用提供了更方便快捷的途徑。第8章MATLAB方程求解8.3.1無約束最優(yōu)化問題求解無約束最優(yōu)化問題的一般描述為：其中，，該數(shù)學(xué)表示的含義亦即求取一組x，使得目標函數(shù)f(x)為最小，故這樣的問題又稱為最小化問題。MATLAB提供了3個求最小值的函數(shù)，它們的調(diào)用格式為：①[x,fval]=fminbnd(filename,x1,x2,option)：求一元函數(shù)在(xl，x2)區(qū)間中的極小值點x和最小值fval。②[x,fval]=fminsearch(filename,x0,option)：基于單純形算法求多元函數(shù)的極小值點x和最小值fval。③[x,fval]=fminunc(filename,x0,option)：基于擬牛頓法求多元函數(shù)的極小值點x和最小值fval。MATLAB沒有專門提供求函數(shù)最大值的函數(shù)，但只要注意到-f(x)在區(qū)間（a，b）上的最小值就是f(x)在（a，b）的最大值，所以fminbnd(-f,x1,x2)返回函數(shù)f(x)在區(qū)間（x1，x2）上的最大值。第8章MATLAB方程求解例8-13求函數(shù)在區(qū)間(-10，-1)和(1，10)上的最小值點。命令如下：>>f=@(x)x-1./x+5;>>[x,fmin]=fminbnd(f,-10,-1)%求函數(shù)在(-10，-1)內(nèi)的最小值點和最小值x=-9.9999fmin=-4.8999>>fminbnd(f,1,10)%求函數(shù)在(1，10)內(nèi)的最小值點ans=1.0001第8章MATLAB方程求解例8-14設(shè)求函數(shù)f在(0.5，0.5，0.5)附近的最小值。建立函數(shù)文件fxyz.m。functionf=fxyz0(u)x=u(1);y=u(2);z=u(3);f=x+y.^2./x/4+z.^2./y+2./z;在MALAB命令行窗口，輸入如下命令：>>[U,fmin]=fminsearch(@fxyz0,[0.5,0.5,0.5])%求函數(shù)的最小值點和最小值U=0.50001.00001.0000fmin=4.0000第8章MATLAB方程求解8.3.2有約束最優(yōu)化問題求解有約束最優(yōu)化問題的一般描述為：其中，，該數(shù)學(xué)表示的含義亦即求取一組x，使得目標函數(shù)f(x)為最小，且滿足約束條件G(x)≤0。記號s.t.是英文subjectto的縮寫，表示x要滿足后面的約束條件。約束條件可以進一步細化為:①線性不等式約束：Ax≤b。②線性等式約束：Aeqx=beq。③非線性不等式約束：C(x)≤0。④非線性等式約束：Ceq(x)=0。⑤x的下界和上界：Lbnd≤x≤Ubnd。第8章MATLAB方程求解MATLAB最優(yōu)化工具箱提供了一個fmincon函數(shù)，專門用于求解各種約束下的最優(yōu)化問題，其調(diào)用格式為：[x,fval]=fmincon(filename,x0,A,b,Aeq,beq,Lbnd,Ubnd,NonF,option)其中，x、fval、filename、x0和option的含義與求最小值函數(shù)相同。其余參數(shù)為約束條件，參數(shù)NonF為非線性約束函數(shù)的M文件名。如果某個約束不存在，則用空矩陣來表示。第8章MATLAB方程求解例8-15求解有約束最優(yōu)化問題。首先編寫目標函數(shù)M文件fop.m。functionf=fop(x)f=0.4*x(2)+x(1)^2+x(2)^2-x(1)*x(2)+1/30*x(1)^3;再設(shè)定約束條件，并調(diào)用fmincon函數(shù)求解此約束最優(yōu)化問題，程序如下：x0=[0.5;0.5];A=[-1,-0.5;-0.5,-1];b=[-0.4;-0.5];lb=[0;0];option=optimset;option.LargeScale='off';option.Display='off';[x,f]=fmincon('fop',x0,A,b,[],[],lb,[],[],option)程序運行結(jié)果為：x=0.33960.3302f=0.2456第8章MATLAB方程求解8.3.3線性規(guī)劃問題求解線性規(guī)劃是研究線性約束條件下線性目標函數(shù)的極值問題的數(shù)學(xué)理論和方法。線性規(guī)劃問題的標準形式為：在MATLAB中求解線性規(guī)劃問題使用函數(shù)linprog，其調(diào)用格式為[x,fval]=linprog(f,A,b,Aeq,beq,Lbnd,Ubnd)其中，x是最優(yōu)解，fval是目標函數(shù)的最優(yōu)值。函數(shù)中的各項參數(shù)是線性規(guī)劃問題標準形式中的對應(yīng)項，x、b、beq、Lbnd、Ubnd是向量，A、Aeq為矩陣，f為目標函數(shù)系數(shù)向量。第8章MATLAB方程求解例8-16求解線性規(guī)劃問題。

命令如下：f=[-5;-4;-6];A=[1,-1,1;3,2,4;3,2,0];b=[20;42;30];Aeq=[];Beq=[];LB=zeros(3,1);[x,favl]=linprog(f,A,b,Aeq,Beq,LB)程序運行結(jié)果為：Optimizationterminated.x=0.000015.00003.0000favl=-78.0000第8章MATLAB方程求解8.4常微分方程初值問題的數(shù)值求解考慮常微分方程的初值問題：y'=f(t,y)，t0≤t≤Ty(t0)=y0所謂其數(shù)值解法，就是求它的解y(t)在節(jié)點t0<t1<…<tm處的近似值y0，y1，…，ym的方法。所求得的y0，y1，…，ym稱為常微分方程初值問題的數(shù)值解。一般采用等距節(jié)點tn=t0+nh，n=0，1，…，m，其中h為相鄰兩個節(jié)點間的距離，叫做步長。常微分方程初值問題的數(shù)值解法多種多樣，比較常用的有歐拉（Euler）法、龍格—庫塔（Runge-Kutta）法、線性多步法、預(yù)報校正法等。第8章MATLAB方程求解8.4.1龍格—庫塔法簡介對于一階常微分方程的初值問題，在求解未知函數(shù)y時，y在t0點的值y(t0)=y0是已知的，并且根據(jù)高等數(shù)學(xué)中的中值定理，應(yīng)有y(t0+h)=y1≈y0+hf(t0,y0)，h>0y(t0+2h)=y2≈y1+hf(t1,y1)一般地，在任意點ti=t0+ih，有：y(t0+ih)=yi≈yi-1+hf(ti-1,yi-1)，i=1，2，…，n當(dāng)（t0，y0）確定后，根據(jù)上述遞推式能計算出未知函數(shù)y在點ti=t0+ih，i=0，1，…，n的一列數(shù)值解：yi=y0，y1，y2，…，yn，i=0，1，…，n當(dāng)然，遞推過程中有一個誤差累計的問題。在實際計算過程中，使用的遞推公式一般進行過改造，著名的龍格—庫塔公式是：y(t0+ih)=yi≈yi-1+(k1+2k2+2k3+k4)其中：k1=f(ti-1,yi-1)k2=f(ti-1+h/2,yi-1+h/2k1)k3=f(ti-1+h/2,yi-1+h/2k2)k4=f(ti-1+h,yi-1+hk3)第8章MATLAB