第五章 數值計算

第五章

數值計算

引言數值計算是MATLAB最強大的功能之一，其矩陣的計算能力強、編程簡單、計算代碼經過精心地優(yōu)化，計算效率很高。在數值計算方面，涉及的數學學科面廣，例如矩陣代數計算、微積分計算、快速傅立葉變換、常微分方程計算、偏微分方程計算、數據擬合等等。MATLAB提供了非常豐富的解決這類問題的命令。

在這一章，我們將主要以實例來介紹解決問題的方法，并假定讀者具有線性代數、數學分析、計算方法等背景知識。具體知識點見導圖5。

5.1矩陣代數的計算

矩陣代數計算包括求矩陣的各種三角分解、線性方程組的求解、矩陣的逆和廣義逆、矩陣的模及條件數、矩陣的譜分解和奇異值分解。矩陣代數的計算是解決數學問題的基礎性方法，任何的數學模型不論是統計模型、微分方程模型還是其他的復雜數學模型最后大部分將進行離散化處理，然后歸為矩陣代數的計算。5.1.1矩陣代數知識點搜索

在參考文獻Reference中進入MATLABFunctionReference，即進入MATLAB的函數參考列表，然后雙擊Mathematics，我們就進入了MATLAB的核心數學方法的函數集（或稱數學方法命令集）。代數計算函數【例5.1】MATLAB在數值計算方面的強大功能。傳統求逆法和MATLAB左除法求線性方程的性能對比。rand('state',12); %選定隨機種子。A=rand(100,100)+1.e8;%生成100×100均勻分布隨機矩陣。x=ones(100,1); %令解向量x為全1的100元列向量。b=A*x;%為使Ax=b方程一致，用A和x生成b向量。cond(A) %求A陣的條件數。ans=1.4426e+012從矩陣的條件數來看，相當大。說明矩陣A是嚴重病態(tài)矩陣，求以它為系數矩陣的線性方程的根將是不精確的，或者說是不穩(wěn)定的。%“求逆”法解恰定方程的誤差、殘差、運算次數和所用時間flops(0);tic%計數器置0，啟動計時器StopwatchTimerxi=inv(A)*b;%xi是用“求逆”法解恰定方程所得的解。ti=toc %關閉計時器，并顯示解方程所用的時間。ci=flops %“求逆”法解方程所用的運算次數eri=norm(x-xi) %解向量xi與真解向量x的范-2誤差。rei=norm(A*xi-b)/norm(b) %方程的范-2相對殘差

ti=0.9300；用時0.93秒ci=2070322；乘除次數eri=3.0708e-004絕對誤差rei=6.6280e-007相對誤差%“左除”法解恰定方程的誤差、殘差、運算次數和所用時間flops(0);tic;xd=A\b; %以上是用“左除”法解方程所得的解。td=toc,cd=flops,erd=norm(x-xd),red=norm(A*xd-b)/norm(b)td=0.2200計算時間cd=741872乘除法的次數erd=3.2243e-004絕對誤差red=2.0095e-016相對誤差對比結果：用左除法是逆陣法運算次數的4/10，時間是逆陣除法的1/4。左除法的相對殘差幾乎為“機器零”，比逆陣法相對精度要高的多。5.1.2矩陣分析MatrixAnalysis矩陣分析主要是考察一個系數矩陣A對線性方程組解的靈敏度，即解的穩(wěn)定性問題。例如對于最小二乘估計的解如果協方差矩陣具有列相關性時，或者說該矩陣接近不滿秩時，我們稱該矩陣為病態(tài)矩陣，它將嚴重影響最小二乘估計的穩(wěn)定性。度量其病態(tài)性的指標為條件數，矩陣的條件數為該矩陣的最大特征值比最小特征值，用來衡量矩陣計算的精確度，當條件數為1時是很好的條件，而條件數非常大時則表明矩陣不滿秩，計算時將產生很大的誤差。矩陣分析包括:條件數的計算、矩陣模的計算、矩陣的行列式的計算和矩陣的標準正交基的計算等等?！纠?.1.1】求矩陣A的條件數、行列式、模等的計算A=[8,-1,-5；-4,4,-2；18,-5,7]%定義一個矩陣c=cond(A)%計算矩陣A的條件數c=8.2546n=norm(A)%計算矩陣的模n=21.5249d=det(A)%計算矩陣的行列式d=412B=orth(A)%求矩陣的一個標準正交基BB=-0.2987-0.9498-0.09330.2469-0.17130.9538-0.92190.26190.2856EYE=B'*B%驗證B是標準正交基EYE=1.0000-0.00000.0000-0.00001.000000.000001.00005.1.3矩陣的分解和三角分解Factorization從計算的角度來看，處理矩陣問題的一個最有效的方法是，將一個一般的矩陣分解為幾個簡單矩陣的乘積形式。上三角矩陣為例，由其性質，兩個上三角陣的和、積仍為上三角陣，上三角陣的特征值就是其對角線元素。系數矩陣為上三角陣的線性方程組只需用回代的方法求解，上三角陣的逆陣仍然是上三角陣。LU分解：將矩陣A分解為A=LU，其中L為單位下三角陣，U為上三角陣。Cholesky分解：將正定對稱陣分解成A=T'T的形式，其中T為正交變換且為上三角陣。QR分解：將矩陣分解為A=QR的形式，其中Q為一個正交陣，R為一個上三角陣?！纠?.1.2】對矩陣A進行LU分解[L,U]=lu(A)%對矩陣A進行LU分解A=8-1-5-44-218-5-7L=0.44440.42311.0000-0.22221.000001.000000U=18.0000-5.0000-7.000002.8889-3.555600-0.3846【例5.1.3】對對稱正定矩陣進行Cholesky分解B=[1-32-310-52-56]T=chol(B)%對對稱矩陣進行Cholesky分解C=T‘*T%對T進行驗證B=1-32-310-52-56T=1-32011001C=1-32-310-52-56

5.1.4求解線性方程組LinearEquations【例5.1.4】求解以下線性方程組

A=[1,-3,2;-3,10,-5;2,-5,6]%定義系數矩陣Ab=[27,-78,64]%定義常數bx=A\b'%求解線性方程組b1=A*x最后的解為：x=(1-47)5.1.5矩陣的特征值或奇異值及其與特征向量EigenvaluesandSingularValues（1）對稱矩陣的譜分解,即求矩陣的特征值和特征向量設矩陣A是N×N維方陣，則譜分解定義如下：其中D為對角矩陣，U為正交矩陣。對角矩陣的元素稱為矩陣A的特征值，矩陣U的列向量為對應于矩陣A的特征值的特征向量。（2）對一般矩陣的奇異值分解設矩陣A為n×m階矩陣，則矩陣的奇異值分解定義為：其中U和V分別為n和m階正交方陣，D為n×m階對角陣，其非對角元素均為0，D的對角線元素稱為矩陣A的奇異值。奇異值的分解與矩陣的譜分解方法類似。求解矩陣A的特征值和特征向量的語法為：[V,D]=eig(A)其中V為矩陣其列向量為矩陣A的特征向量，D為對角矩陣對角線元素為矩陣A的特征值?！纠?.1.5】求矩陣的特征值和奇異值。（1）求上面矩陣A的特征值和特征向量。[V,D]=eig(A)%求矩陣A的特征值和特征向量V=0.9654-0.00910.26060.22110.5581-0.7998-0.13810.82970.5408D=0.02660002.615100014.3583（2）求下列矩陣B的奇異值及其對應的向量。B=[94；68；27]%矩陣B為長方形矩陣[U,S,V]=svd(B)%求矩陣B的奇異值和對應的兩個向量B=946827U=-0.61050.71740.3355-0.6646-0.2336-0.7098-0.4308-0.65630.6194S=14.9359005.188300V=-0.69250.7214-0.7214-0.69255.2多項式與插值在數學建模競賽中或在科學研究中，使用頻率最高的數據加工技術是對實驗數據的擬合，擬合有很多方法如多項式方法、最小二乘法和樣條插法。當我們對給定的數據利用上面的方法選擇了擬合方法后，對任意給定的自變量點（插值點）我們可以計算該點的函數值了。

MATLAB約定多項式如下：我們可以用行向量來表示一個多項式。如，。我們可以進行多項式的乘、除、因子的提取、多項式的簡化、求多項式的根等等常規(guī)的多項式的計算和多項式的擬合。注意：如果某系數為0，必須寫上以占一個位置。5.2.1多項式的創(chuàng)建方法（1）直接用行向量創(chuàng)建，如：P=[1，2，3，5，6，7]%這創(chuàng)建了一個5階多項式。

（2）利用命令：P=ploy（A）%產生多項式系數向量。這里A為方陣，即多項式P為A的特征多項式?！纠?.2.1】求三階魔方方陣M的特征多項式。M=magic(3)PM=poly(M)； %A的特征多項式PPM=poly2str(PM,'x') %將多項式轉成習慣的方式RM=roots(PM)%求A或者說特征多項式的根PPA=x^3-15x^2-24x+360RM=15.0000-4.89904.89905.2.2計算多項式的操作函數與例題多項式的運算包括，多項式相乘、多項式相除、求多項式的公因式等等運算，表5.2.1列出一些基本多項式運算命令【例5.2.2】由給定的根向量求多項式的系數向量。R=[-0.5,-0.3+0.4*i,-0.3-0.4*i]; %根向量P=poly(R) %R的特征多項式PR=real(P) %求PR的實部PPR=poly2str(PR,'x')計算結果如下：P=1.00001.10000.55000.1250PR=1.00001.10000.55000.1250PPR=x^3+1.1x^2+0.55x+0.125注意：（1）要形成實系數多項式，根向量中的復數必須共軛成對。（2）含復數的根向量產生的多項式系數向量P有可能存在截斷誤差級的虛部。我們可用real命令把虛數過濾掉。【例5.2.3】計算有理式的“商”及“余”多項式p1=conv([1,0,2],conv([1,4],[1,1])); %計算分子多項式p2=[1011]; %注意缺項補零[q,r]=deconv(p1,p2);cq='商多項式為';cr='余多項式為';disp([cq,poly2str(q,'s')])disp([cr,poly2str(r,'s')])計算結果為：商多項式為s+5余多項式為5s^2+4s+3【例5.2.4】我們對中國1952年到1998年的國民總收入GDP的數據，用1、2、3和4階多項式來擬合，并將結果用圖形表示出來。GDP的數據文件為gdp.txtclear;AA=load('e:\data\gdp.txt');x=AA(:,1);y=AA(:,2);p1=polyfit(x,y,1)%用一階多項式擬合p2=polyfit(x,y,2)%用二階多項式擬合p3=polyfit(x,y,4)%用四階多項式擬合p4=polyfit(x,y,6)%用六階多項式擬合f1=polyval(p1,x);%計算一階多項式的插值點f2=polyval(p2,x);%計算二階多項式的插值點f3=polyval(p4,x);%計算三階多項式的插值點f4=polyval(p4,x);%計算四階多項式的插值點subplot(2,2,1),plot(x,y,'.',x,f1,'-'),title('一階圖')subplot(2,2,2),plot(x,y,'.',x,f2,'-'),title('二階圖')subplot(2,2,3),plot(x,y,'.',x,f3,'-'),title('三階圖')subplot(2,2,4),plot(x,y,'.',x,f4,'-'),title('四階圖')從結果圖上看，四階多項式和六階多項式擬合的好。我們可以取四階多項式，要注意的是，并不是階數越多越好。5.2.4插值計算插值與多項式擬合的不同之處在于，對給定的點多項式擬合不必點點通過，而插值建立的函數是點點通過給定的數據點。MATLAB提供了豐富的插值方法，并配有樣條插值工具箱。5.2.5一維插值多項式一維插值多項式的基本語法為：yi=interp1(x,y,xi,yi,method)這里x，y為給定的數據，利用該數據建立一維插值多項式，xi為待計算的插值點，yi為插值點的函數值。method：為建立插值多項式的方法，可取method='nearest'：利用最近臨方法建立插值多項式method='linear'：利用線性方法建立插值多項式method='spline'：利用樣條函數方法建立插值多項式method='pchip'：三階厄米特方法建立逐段樣條插值多項式method='cubic'：利用三階樣條函數建立插值多項式【例5.2.5】對國民經濟數據從1955年以5年為間隔建立插值函數，并在求1955年、1986年和預測1996年的數據。B=ones(3,2)A=load('e:\data\zggdp.txt')x=A(:,1)y=A(:,2)ly=log(y)x1=x([4,9,14,19,24,29,34,39,44])y1=y([4,9,14,19,24,29,34,39,44])xx=[1955,1986,1996]yy=interp1(x,y,xx,'spline')B(:,1)=xx';B(:,2)=yy';B結果為：1955910196610201199666851【例5.2.6】在極坐標下插值的例子，給定橢圓上的四個點，用樣條插值并作圖。theta=[0:0.5:2]*pi; %產生自變量四個樣點y=[-0.51-0.5-10.51-0.50.510.5-1-0.510.5]; %給出相應的因變量的點和初值theta2=linspace(theta(1),theta(end),100);%參量稠密化yy=spline(theta,y,theta2); %求稠密點上的插值plot(yy(1,:),yy(2,:),'b',y(1,:),y(2,:),'or')5.2.6高維插值我們以二維插值為例，其語法為：ZI=interp2(X,Y,Z,XI,YI,method)其中的參數和一維的類似?！纠?.2.7】以peaks函數為例，首先取較少的數據。然后對其求插值函數，并比較各方法的差異。[x,y]=meshgrid(-3:1:3);z=peaks(x,y);surf(x,y,z);title('較少點的表面圖')[xi,yi]=meshgrid(-3:0.25:3);zi1=interp2(x,y,z,xi,yi,'nearest');subplot(2,2,1),surf(xi,yi,zi1),title('nearest方法')zi2=interp2(x,y,z,xi,yi,'bilinear');subplot(2,2,2),surf(xi,yi,zi2),title('bilinear方法')zi3=interp2(x,y,z,xi,yi,'cubic');subplot(2,2,3),surf(xi,yi,zi3),title('bicubict方法')zi4=interp2(x,y,z,xi,yi);subplot(2,2,4),surf(xi,yi,zi4),title('系統內定方法')5.3微積分在大學生數學建模競賽中微積分知識出現的比率是很高的，也是很基礎知識。比如2002年大學生數模競賽中，A題就是一個有關微積分知識的應用題，即《車燈線光源的優(yōu)化設計》的問題。我們知道車燈是一個拋物線型結構，光源在旋轉拋物面的焦點處且形狀為小長條形。我們要研究光源條的長短對燈光照射的影響等等。5.3.1函數的數值導數與切平面函數曲面的法線用i，j，k表示直角坐標系的單位向量，曲面在（x0,y0）處的法線（Normal）n是

在MATLAB中計算與繪制法線的命令為：[NX，NY，NZ]=surfnorm（X，Y，Z）在MATLAB中計算與繪制法線的命令為：[NX，NY，NZ]=surfnorm（X，Y，Z）【例5.3.1】曲面法線演示。y=-1:0.1:1;x=2*cos(asin(y));%旋轉曲面的“母線”[X,Y,Z]=cylinder(x,20); %形成旋轉曲面surfnorm(X(:,11:21),Y(:,11:21),Z(:,11:21)); %在曲面上畫法線view([120,18]) %控制觀察角5.3.2偏導數和梯度函數的偏導數定義如下：全微分可借助梯度（Gradient）定義為：數值差分命令的基本語法：DX=diff（X）求X相鄰元素的一階差分DX=diff（X，n）求X相鄰元素的n階差分DX=diff（X，n，dim）命令維dim上,求相鄰元素的n階差分差分的含義為：設序列則一階差分的定義為：n階差分依次類推?！纠?.3.2】差分方法應用十分廣泛，例如在經濟數據中一階差分表示發(fā)展速度，二階差分表示經濟發(fā)展的加速度。對我國52~98年國民經濟數據GDP發(fā)展趨勢進行分析，求其一、二階差分。并作圖。A=load('e:\data\gdp.txt')B=A(:,2)B1=diff(B);B2=diff(B,2)subplot(1,3,1),plot(B),title('發(fā)展趨勢')subplot(1,3,2),plot(B1),title('發(fā)展速度')subplot(1,3,3),plot(B2),title('發(fā)展加速度')梯度:基本語法：[FX，FY]=gradient（F，h）求二元函數的梯度[FX，FY，FZ，…]=gradient（F，h1,h2,…）【例5.3.3】對函數作梯度圖v=-2:0.2:2;[x,y]=meshgrid(v);z=x.*exp(-x.^2-y.^2);subplot(1,2,1);surf(x,y,z),subplot(1,2,2);[px,py]=gradient(z,.2,.2);contour(x,y,z),holdon,quiver(x,y,px,py)方向法線的繪圖命令：基本語法為：quiver3(Z,U,V,W)quiver3(X,Y,Z,U,V,W)quiver3(...,scale)quiver3(...,LineSpec)這里：X，Y，Z：為曲線上的點。U，V，W：為每一點上切面的法線方向。scale：為實數值，可正可負。表示方向箭頭的指向?！纠?.3.4】求函數的曲面圖，且作出每點的切面法向量圖形。[X,Y]=meshgrid(-2:0.25:2,-1:0.2:1);%產生網格點Z=X.*exp(-X.^2-Y.^2);%計算每個網格點的函數值[U,V,W]=surfnorm(X,Y,Z);%計算每個點的法向量分量quiver3(X,Y,Z,U,V,W,0.5);%作法向量的圖形holdonsurf(X,Y,Z);colormaphsvview(-35,45)axis([-22-11-.6.6])holdoff【例5.3.5】綜合應用。我們來作一個旋轉拋物面上某點的切面和該切面的法線。這里用到了作方向法線的命令。x=0:0.1:1;y=sqrt(2*x);%旋轉曲面的“母線”[X,Y,Z]=cylinder(y,20); %形成旋轉曲面[U,V,W]=surfnorm(X,Y,Z);Z1=Z(4)-(U(4)*(X-X(4))+V(4)*(Y-Y(4)))/W(4)%求點X(4),Y(4),Z(4)的切平面surf(X,Y,Z); %在曲面上畫法線holdon;surf(X,Y,Z1)；holdon在上面點處畫切平面plot3(X(4),Y(4),Z(4),'r.','MarkerSize',30)%將該點突出放大并用紅色表示U(1:3)=NaN;U(5:end)=NaNholdon;quiver3(X,Y,Z,U,V,W,-4)%作該點處的法線，長度為4,反方向axisequal%該命令是必須的，使得視覺正確5.3.3數值積分對于形如的積分積分的一般命令為：q=quad(fun,a,b)q=quad(fun,a,b,tol)q=quad(fun,a,b,tol,trace)q=quad(fun,a,b,tol,trace,p1,p2,...)[q,fcnt]=quadl(fun,a,b,...)其中fun表示被積函數，系數a，b為積分的上下限，tol為用戶給定的誤差限，如不給出則系統內定的值為0.0000001?！纠?.3.6】以下是一個比較著名的例子，計算積分Q=quad('1./(1+25*x.^2)',-1,1)

ans=0.5494也可以使用內聯函數定義被積函數F=inline('1./(1+25*x.^2)');%定義一個內聯函數Q=quad(F,-1，1);ans=0.5494【例5.3.7】對空間曲線求積分，設空間曲線的參數方程為：x=t.*sin(t),y=t.*cos(t),z=tt=0:pi/50:10*pi;plot3(t.*sin(t),t.*cos(t),t,'LineWidth',2)f=inline('sqrt(4*(t.*sin(t)).^2+(t.*cos(t)).^2…+1)');len=quad(f,0,pi)計算結果為：len=8.51915.3.4二重積分二重積分的原理與一重積分一樣，使用也很方便。其語法為：q=dblquad(fun,xmin,xmax,ymin,ymax)q=dblquad(fun,xmin,xmax,ymin,ymax,tol)q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)【例3.3.8】求積分q=dblquad('exp(-x.^2.-y^2)',-1,1,-1,1)

q=2.23105.4求函數的根函數求根或稱求函數的根，是初等數學就學習過的。然而對于復雜的函數、方程組求根的運算是用搜索的方法根據函數的特性逐步接近某個根，在使用MATLAB命令時應該對函數有所了解，給出某個根附近的初始點或根所在的某一個區(qū)間。主要求根命令為：roots：求多項式的根fzero：求一般一元函數的根fsolve：求非線性方程組的根5.4.1求多項式和一般一元函數的根【例5.4.1】求多項式的根clear;p=[1,-6,-72,-27]r=roots(p)p=1-6-72-27r=12.1229

-5.7345

-0.3884求一元函數的0點命令fzero其基本語法為：x=fzero(fun,x0)x=fzero(fun,x0,options)x=fzero(fun,x0,options,P1,P2,...)[x,fval]=fzero(...)[x,fval,exitflag]=fzero(...)[x,fval,exitflag,output]=fzero(...)這里fun：被求解的函數，可用inline或字符串定義x0：為初初始值，可以是一個數，也可以是一個區(qū)間fval：目標函數值options：輸出顯示項，可選擇display：'off':沒有輸出'iter'：將每一步的迭代結果都輸出來'notify'：只有當迭代不收斂時才輸出警告信息（系統內定值）TolX：迭代終止精度。【例5.4.2】在[1，2]范圍sinx^2的根并顯示詳細迭代步驟。options=optimset('Display','iter')x=fzero('sin(x*x)',[12],options)結果為：Func-countxf(x)Procedure110.841471initial22-0.756802initial31.526490.725271interpolation41.758210.0502804interpolation51.77439-0.00687537interpolation61.772453.0226e-005interpolation71.772451.62826e-008interpolation81.77245-1.2098e-015interpolation91.772451.89882e-015interpolationZerofoundintheinterval:[1,2].x=1.7725對于任意f(x)可能有0點，可能沒有，可能有多個0點，甚至可能有無窮的0點。因此找不到一個通解，我們可以通過人機對話，將函數顯示在屏幕上，在可能的點上用鼠標采樣得到初始點，然后讓計算機根據這些點自動分別地求出每一個根?！纠?.4.3】編制程序，對任何函數求根，在運行該程序時將顯示函數與零點的圖形，用戶在可能的零點附近用鼠標點擊采點，程序通過循環(huán)計算每個采集點的根，并輸出結果。以函數

求零點為例，綜合敘述相關命令的用法。y=inline('sin(t)^2*exp(-a*t)-b*abs(t)','t','a','b'); %建立內聯函數a=0.1;b=0.5;t=-10:0.01:10;%對自變量采樣，采樣步長不宜太大。y_char=vectorize(y);Y=feval(y_char,t,a,b); clf%作人機對話界面plot(t,Y,'r');holdon,plot(t,zeros(size(t)),'k'); %畫坐標橫軸xlabel('t');ylabel('y(t)'),holdofftitle('用鼠標在函數零點附近連續(xù)點五點')%循環(huán)使用求解命令fzero，求鼠標采集點附近的函數根。fori=1:5[t1,y1]=ginput(1);[ttt(i),yyy(i),exitflag]=fzero(y,t1,[],a,b) %在ttt（i）附近搜索0點endA(:,1)=ttt';A(:,2)=yyy';A則計算結果為：A=xf（x）-2.00740.0000-2.00740.0000-0.51980.00001.67380.00001.6738-0.00005.4.2求多元非線性方程組的根求多元非線性方程組的根，算法較為復雜，我們可以利用fsolve命令來求解，其命令語法為：x=fsolve(fun,x0)x=fsolve(fun,x0,options)x=fsolve(fun,x0,options,P1,P2,...)[x,fval]=fsolve(...)[x,fval,exitflag]=fsolve(...)[x,fval,exitflag,output]=fsolve(...)[x,fval,exitflag,output,jacobian]=fsolve(...)這里fun：多元非線性方程組x0：一向量，初始值fval：目標函數值output:輸出內容控制iterations：迭代次數輸出funcCount：計算函數值algorithm：算法使用cgiterations：PCG迭代次數（僅對大型問題）stepsize：：最終步長的選擇（僅對中型問題）firstorderoptoptions：輸出顯示項，可選擇Diagnostics：輸出診斷信息display：'off':沒有輸出'iter'：將每一步的迭代結果都輸出來'notify'：只有當迭代不收斂時才輸出警告信息（系統內定值）Jacobian:雅克比方法'on':精確雅克比方法'off'：近似雅克比方法MaxFunEvals：函數的最大上限MaxIter：最大迭代次數TolFun：終止容忍函數值TolX：終止迭代容忍度exitflag：退出條件>0表示函數收斂到某解，正常退出達到最大給定迭代次數，退出<0求解過程不收斂【例5.4.4】求解二元函數方程組

在[-5，-5]為初始值的零點。F='[2*x(1)-x(2)-exp(-x(1)),-x(1)+2*x(2)-exp(-x(2))]';x0=[-5;-5];%初始值options=optimset('Display','iter');%為Option設置值[x,fval]=fsolve(F,x0,options)%求零點x=0.56710.5671fval=1.0e-008*-0.5319-0.53195.5求函數極值及最優(yōu)化問題MATLAB提供一些求函數極值的函數，如單變量函數求極小fmins和多變量函數求極小fminsearch等。這兩個函數的語法和參數與求根函數類似不再綴述。我們也可以打開最優(yōu)化工具箱使用里面的更專業(yè)的方法。現將主要的命令列表如下：5.5.1無約束最優(yōu)化問題【例5.5.1】求的最小值，給定初始值-1[x,fval,exitflag]=fminbnd('x^2',-1,'x',…optimset('TolX',1e-12,'Display','off'))結果為：x=0fval=0exitflag=1【例5.5.2】一個典型的函數是Rosenbrockbanana函數，其典型的最小值坐標我們已經知道為（1，1）。首先將它的圖形作出。（1）作Rosenbrockbanana函數的三維等高圖x=-3:0.1:3;y=-2:0.1:4;%[-3，3]×[-2，4]區(qū)域取點[X,Y]=meshgrid(x,y);%根據點在區(qū)域打上網格F=100*(Y-X.^2).^2+(1-X).^2;%每一點處計算函數值contour3(X,Y,F,300),%作三維等高線圖形xlabel('x'),ylabel('y'),axis([-3,3,-2,4,0,inf]),view([161,45])holdon,plot3(1,1,0,'.r','MarkerSize',20),holdoff我們可以看到有一條香蕉狀的谷底，事實上該函數的最小值不是對所有的方法都能找到的，因為香蕉狀的谷底較平坦，因此Rosenbrockbanana函數成了著名的搜索方法的測試函數。（2）建立用于搜索的Rosenbrockbanana函數的符號形式ff=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2','x');

（3）用單純形法求極小值點x0=[-1.2,1];

%給出X-Y片面的初始搜索點[sx,sfval,sexit,soutput]=fminsearch(ff,x0)

%進行搜索Optimizationterminatedsuccessfully:thecurrentxsatisfiestheterminationcriteriausingOPTIONS.TolXof1.000000e-004andF(X)satisfiestheconvergencecriteriausingOPTIONS.TolFunof1.000000e-004最優(yōu)點為：sx=1.00001.0000最優(yōu)值為：sfval=8.1777e-010sexit=1soutput=iterations:85funcCount:159algorithm:'Nelder-Meadsimplexdirectsearch'5.5.2有約束最優(yōu)化問題(1)線性規(guī)劃問題模型為：（5.5.1）線性規(guī)劃問題命令的一般語法為：x=linprog(f,A,b,Aeq,beq)x=linprog(f,A,b,Aeq,beq,lb,ub)x=linprog(f,A,b,Aeq,beq,lb,ub,x0)x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options)這里：f：目標函數系數向量A：不等式約束矩陣b：不等式約束右邊項系數向量Aeq：等式約束矩陣beq：等式約束右邊項系數向量lb：自變量下限向量ub：自變量上限向量[x,fval]=linprog(...)[x,fval,exitflag]=linprog(...)[x,fval,exitflag,output]=linprog(...)[x,fval,exitflag,output,lambda]=linprog(...)【例5.5.3】求以下線性規(guī)劃問題的解（5.5.2）編程如下：f=[-5;-4;-6]A=[1-11；324；320];b=[20;42;30];lb=zeros(3,1);[x,fval,exitflag,output,lambda]=linprog(f,A,b,[],[],lb)f=-5-4-6Optimizationterminatedsuccessfully.x=0.000015.00003.0000

最優(yōu)解fval=-78.0000exitflag=1output=iterations:6cgiterations:0algorithm:'lipsol'lambda=ineqlin:[3x1double]eqlin:[0x1double]upper:[3x1double]lower:[3x1double](1)二次規(guī)劃二次規(guī)劃模型與線性規(guī)劃類似，只是其目標函數具有變量的二次項。其模型為：(5.5.3)二次規(guī)劃命令語法為：x=quadprog(H,f,A,b)x=quadprog(H,f,A,b,Aeq,beq)x=quadprog(H,f,A,b,Aeq,beq,lb,ub)x=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0)x=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options)x=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options,p1,p2,...)[x,fval]=quadprog(...)[x,fval,exitflag]=quadprog(...)[x,fval,exitflag,output]=quadprog(...)[x,fval,exitflag,output,lambda]=quadprog(...)這里H為變量二次型矩陣，其他參數同線性規(guī)劃類似【例5.5.4】求以下二次規(guī)劃的最優(yōu)值與最優(yōu)解。目標函數其中%首先數據準備H=[1-1;-12];f=[-2;-6]A=[11;-12;21];b=[2;2;3]lb=zeros(2,1)[x,fval,exitflag,output,lambda]=quadprog(H,f,A,b,[],[],lb)%求二次規(guī)劃x0=[10;10;10];%StartingguessatthesolutionA=[-1,-2,-2;1,2,2]b=[0;72][x,fval]=fmincon(f,x0,A,b)計算結果為：Optimizationterminatedsuccessfully.x=0.66671.3333最優(yōu)解fval=-8.2222最優(yōu)值exitflag=1output=iterations:3algorithm:'medium-scale:active-set'firstorderopt:[]cgiterations:[]lambda=lower:[2x1double]upper:[2x1double]eqlin:[0x1double]ineqlin:[3x1double]5.5.3非線性有約束最優(yōu)化問題對于非線性函數的求最優(yōu)化問題情況較為復雜，沒有一個固定的方法。非線性函數的求最優(yōu)化模型為：其中f(x)，c(x)，ceq(x)，為非線形函數。求非線性函數的求最優(yōu)化問題命令的語法為：x=fmincon(fun,x0,A,b)x=fmincon(fun,x0,A,b,Aeq,beq)x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub)x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2,...)[x,fval]=fmincon(...)[x,fval,exitflag]=fmincon(...)[x,fval,exitflag,output]=fmincon(...)[x,fval,exitflag,output,lambda]=fmincon(...)[x,fval,exitflag,output,lambda,grad]=fmincon(...)[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(...)這里fun為非線性函數?！纠?.5.5】求以下非線性函數的最小值及最優(yōu)解。搜索初始值為[10，10，10]'解：這里x0=[10;10;10];%StartingguessatthesolutionA=[-1,-2,-2;1,2,2]b=[0;72][x,fval]=fmincon(f,x0,A,b)計算結果為：x=24.0012.0012.0fval=-34565.5.4多目標線性規(guī)劃在實際問題中目標函數可能有多個，例如一個大型的工程其最終的目標是多方面的，既要考慮工程的成本盡量少，又要考慮企業(yè)投產帶來的生態(tài)和社會效應等多方面的因素。這就是一個多目標的規(guī)劃問題。多目標規(guī)劃命令的語法為：x=fgoalattain(fun,x0,goal,weight)x=fgoalattain(fun,x0,goal,weight,A,b)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon,options)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2,...)[x,fval]=fgoalattain(...)[x,fval,attainfactor]=fgoalattain(...)[x,fval,attainfactor,exitflag]=fgoalattain(...)其中goal：為各目標函數要達到的上界，goal為向量，維數是目標的個數。weight：為搜索步長權重?！纠?.5.6】求以下多目標規(guī)劃。我們編程序如下：f='[x(1)+2*x(2);-x(1)-x(2)]'goal=[2,-1];weight=[1,1];x0=[2,1];A=[1,0;0,1];b=[1;1]lb=zeros(2,1);[x,fval,attainfactor,exitflag]=…fgoalattain(f,x0,goal,weight,A,b,[],[],lb,[])計算結果為：最優(yōu)解為：x=1.00000.3333最終多目標值：fval=1.6667-1.33335.6傅立葉分析與快速傅立葉變換傅立葉分析在連續(xù)和離散脈沖信號處理上是非常有用的工具，它可以將輸入的波長進行變換，分離出各種不同的頻率或者提取輸入波長不同的數字特征。一般情況下，我們都是對連續(xù)的波進行離散化處理，得到一個離散的向量進行傅立葉變換，簡記為（DFT）。而快速傅立葉變換（FTT）是處理離散傅立葉變換最強有力的工具，非常適合處理諸如時間序列、信號和圖象序列數據的分析，可以對信號進行濾波處理、變換、譜分析和估計等等?？梢詫⒏盗⑷~分析看成是一個過濾器。例如：

x1,x2,…,xN

X1,X2,…,XN

X1,X2,…,XN

x1,x2,…,Xn傅立葉變換公式為：DFTIDFT傅立葉的逆變換公式為：如果x（n）為實數，則可寫成這里【例5.6.1】這是時間序列課程中一個著名的例子，觀測太陽黑子變化規(guī)律。我們已經有300多年有關太陽黑子的數據，稱為著名的Wolfernumber。并且我們也已經知道太陽黑子的變化周期是11年，本例用fft方法對該時間序列進行變換，并進行周期特征的提取。loadsunspot.dat%讀入太陽黑子數據year=sunspot(:,1);wolfer=sunspot(:,2);plot(year,wolfer)title('300年的太陽黑子數據圖')我們可以看出太陽黑子是呈周期性變化的，我們初步觀察它的周期是在10~12年之間。為了得出更精確的周期，我們進行序列的特征提取，即進行fft變換。Y=fft(wolfer)%對序列進行快速傅立葉變換通過變換Y的元素已為復數，我們來尋找最大周期。首先求Y元素的模長平方，記為power，建立周期和power的圖形，然后找出所求的周期。N=length(Y);Y(1)=[]power=abs(Y(1:N/2)).^2;nyquist=1/2;freq=(1:N/2)/(N/2)*nyquist;period=1./freq;plot(period,power),gridon,axis([04002e7])ylabel('Power')xlabel('Period(Years/Cycle)')title('周期提取示意圖')從圖中我們看出power最高點處周期的值大約為11年。為精確起見，我們來搜索最高點處的周期值。[mp,index]=max(power);period(index)ans=11.0769即周期為11.0769年5.7求解常微分方程

常微分方程是這樣一類數學問題，即其中某個變量和其他變量之間的函數依賴關系是未知的，但是這個未知的函數關系以及它的某些階的導數（或函數）連同自變量都由一個已知的方程聯系在一起，這樣的方程稱為微分方程。如果未知函數是一元的，那么對應的微分方程稱為常微分方程（OrdinaryDifferentialEquations，ODEs）；如果未知函數是多元的，那么對應的微分方程稱為偏微分方程。常微分方程的求解問題分為：初值問題（InitialValueProblem，VIP），和邊值問題（BoundaryValueProblem，BVP），以及求解微分代數方程組問題（differential-algebraicequationsDAEs）。相比而言邊值問題難度更大，且種類繁多。只有很少一類常微分問題可以得到解析解，這里介紹的方法都是常微分方程的數值解問題，數值解沒有一個統一的方法，要視問題的具體情況具體解決。初值問