第五章 數(shù)值計(jì)算

第五章

數(shù)值計(jì)算

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

在這一章，我們將主要以實(shí)例來介紹解決問題的方法，并假定讀者具有線性代數(shù)、數(shù)學(xué)分析、計(jì)算方法等背景知識。具體知識點(diǎn)見導(dǎo)圖5。

5.1矩陣代數(shù)的計(jì)算

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

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

ti=0.9300；用時(shí)0.93秒ci=2070322；乘除次數(shù)eri=3.0708e-004絕對誤差rei=6.6280e-007相對誤差%“左除”法解恰定方程的誤差、殘差、運(yùn)算次數(shù)和所用時(shí)間flops(0);tic;xd=A\b; %以上是用“左除”法解方程所得的解。td=toc,cd=flops,erd=norm(x-xd),red=norm(A*xd-b)/norm(b)td=0.2200計(jì)算時(shí)間cd=741872乘除法的次數(shù)erd=3.2243e-004絕對誤差red=2.0095e-016相對誤差對比結(jié)果：用左除法是逆陣法運(yùn)算次數(shù)的4/10，時(shí)間是逆陣除法的1/4。左除法的相對殘差幾乎為“機(jī)器零”，比逆陣法相對精度要高的多。5.1.2矩陣分析MatrixAnalysis矩陣分析主要是考察一個系數(shù)矩陣A對線性方程組解的靈敏度，即解的穩(wěn)定性問題。例如對于最小二乘估計(jì)的解如果協(xié)方差矩陣具有列相關(guān)性時(shí)，或者說該矩陣接近不滿秩時(shí)，我們稱該矩陣為病態(tài)矩陣，它將嚴(yán)重影響最小二乘估計(jì)的穩(wěn)定性。度量其病態(tài)性的指標(biāo)為條件數(shù)，矩陣的條件數(shù)為該矩陣的最大特征值比最小特征值，用來衡量矩陣計(jì)算的精確度，當(dāng)條件數(shù)為1時(shí)是很好的條件，而條件數(shù)非常大時(shí)則表明矩陣不滿秩，計(jì)算時(shí)將產(chǎn)生很大的誤差。矩陣分析包括:條件數(shù)的計(jì)算、矩陣模的計(jì)算、矩陣的行列式的計(jì)算和矩陣的標(biāo)準(zhǔn)正交基的計(jì)算等等?！纠?.1.1】求矩陣A的條件數(shù)、行列式、模等的計(jì)算A=[8,-1,-5；-4,4,-2；18,-5,7]%定義一個矩陣c=cond(A)%計(jì)算矩陣A的條件數(shù)c=8.2546n=norm(A)%計(jì)算矩陣的模n=21.5249d=det(A)%計(jì)算矩陣的行列式d=412B=orth(A)%求矩陣的一個標(biāo)準(zhǔn)正交基BB=-0.2987-0.9498-0.09330.2469-0.17130.9538-0.92190.26190.2856EYE=B'*B%驗(yàn)證B是標(biāo)準(zhǔn)正交基EYE=1.0000-0.00000.0000-0.00001.000000.000001.00005.1.3矩陣的分解和三角分解Factorization從計(jì)算的角度來看，處理矩陣問題的一個最有效的方法是，將一個一般的矩陣分解為幾個簡單矩陣的乘積形式。上三角矩陣為例，由其性質(zhì)，兩個上三角陣的和、積仍為上三角陣，上三角陣的特征值就是其對角線元素。系數(shù)矩陣為上三角陣的線性方程組只需用回代的方法求解，上三角陣的逆陣仍然是上三角陣。LU分解：將矩陣A分解為A=LU，其中L為單位下三角陣，U為上三角陣。Cholesky分解：將正定對稱陣分解成A=T'T的形式，其中T為正交變換且為上三角陣。QR分解：將矩陣分解為A=QR的形式，其中Q為一個正交陣，R為一個上三角陣。【例5.1.2】對矩陣A進(jìn)行LU分解[L,U]=lu(A)%對矩陣A進(jìn)行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】對對稱正定矩陣進(jìn)行Cholesky分解B=[1-32-310-52-56]T=chol(B)%對對稱矩陣進(jìn)行Cholesky分解C=T‘*T%對T進(jìn)行驗(yàn)證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]%定義系數(shù)矩陣Ab=[27,-78,64]%定義常數(shù)bx=A\b'%求解線性方程組b1=A*x最后的解為：x=(1-47)5.1.5矩陣的特征值或奇異值及其與特征向量EigenvaluesandSingularValues（1）對稱矩陣的譜分解,即求矩陣的特征值和特征向量設(shè)矩陣A是N×N維方陣，則譜分解定義如下：其中D為對角矩陣，U為正交矩陣。對角矩陣的元素稱為矩陣A的特征值，矩陣U的列向量為對應(yīng)于矩陣A的特征值的特征向量。（2）對一般矩陣的奇異值分解設(shè)矩陣A為n×m階矩陣，則矩陣的奇異值分解定義為：其中U和V分別為n和m階正交方陣，D為n×m階對角陣，其非對角元素均為0，D的對角線元素稱為矩陣A的奇異值。奇異值的分解與矩陣的譜分解方法類似。求解矩陣A的特征值和特征向量的語法為：[V,D]=eig(A)其中V為矩陣其列向量為矩陣A的特征向量，D為對角矩陣對角線元素為矩陣A的特征值。【例5.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的奇異值及其對應(yīng)的向量。B=[94；68；27]%矩陣B為長方形矩陣[U,S,V]=svd(B)%求矩陣B的奇異值和對應(yīng)的兩個向量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多項(xiàng)式與插值在數(shù)學(xué)建模競賽中或在科學(xué)研究中，使用頻率最高的數(shù)據(jù)加工技術(shù)是對實(shí)驗(yàn)數(shù)據(jù)的擬合，擬合有很多方法如多項(xiàng)式方法、最小二乘法和樣條插法。當(dāng)我們對給定的數(shù)據(jù)利用上面的方法選擇了擬合方法后，對任意給定的自變量點(diǎn)（插值點(diǎn)）我們可以計(jì)算該點(diǎn)的函數(shù)值了。

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

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

在MATLAB中計(jì)算與繪制法線的命令為：[NX，NY，NZ]=surfnorm（X，Y，Z）在MATLAB中計(jì)算與繪制法線的命令為：[NX，NY，NZ]=surfnorm（X，Y，Z）【例5.3.1】曲面法線演示。y=-1:0.1:1;x=2*cos(asin(y));%旋轉(zhuǎn)曲面的“母線”[X,Y,Z]=cylinder(x,20); %形成旋轉(zhuǎn)曲面surfnorm(X(:,11:21),Y(:,11:21),Z(:,11:21)); %在曲面上畫法線view([120,18]) %控制觀察角5.3.2偏導(dǎo)數(shù)和梯度函數(shù)的偏導(dǎo)數(shù)定義如下：全微分可借助梯度（Gradient）定義為：數(shù)值差分命令的基本語法：DX=diff（X）求X相鄰元素的一階差分DX=diff（X，n）求X相鄰元素的n階差分DX=diff（X，n，dim）命令維dim上,求相鄰元素的n階差分差分的含義為：設(shè)序列則一階差分的定義為：n階差分依次類推?！纠?.3.2】差分方法應(yīng)用十分廣泛，例如在經(jīng)濟(jì)數(shù)據(jù)中一階差分表示發(fā)展速度，二階差分表示經(jīng)濟(jì)發(fā)展的加速度。對我國52~98年國民經(jīng)濟(jì)數(shù)據(jù)GDP發(fā)展趨勢進(jìn)行分析，求其一、二階差分。并作圖。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，F(xiàn)Y]=gradient（F，h）求二元函數(shù)的梯度[FX，F(xiàn)Y，F(xiàn)Z，…]=gradient（F，h1,h2,…）【例5.3.3】對函數(shù)作梯度圖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：為曲線上的點(diǎn)。U，V，W：為每一點(diǎn)上切面的法線方向。scale：為實(shí)數(shù)值，可正可負(fù)。表示方向箭頭的指向。【例5.3.4】求函數(shù)的曲面圖，且作出每點(diǎn)的切面法向量圖形。[X,Y]=meshgrid(-2:0.25:2,-1:0.2:1);%產(chǎn)生網(wǎng)格點(diǎn)Z=X.*exp(-X.^2-Y.^2);%計(jì)算每個網(wǎng)格點(diǎn)的函數(shù)值[U,V,W]=surfnorm(X,Y,Z);%計(jì)算每個點(diǎn)的法向量分量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】綜合應(yīng)用。我們來作一個旋轉(zhuǎn)拋物面上某點(diǎn)的切面和該切面的法線。這里用到了作方向法線的命令。x=0:0.1:1;y=sqrt(2*x);%旋轉(zhuǎn)曲面的“母線”[X,Y,Z]=cylinder(y,20); %形成旋轉(zhuǎn)曲面[U,V,W]=surfnorm(X,Y,Z);Z1=Z(4)-(U(4)*(X-X(4))+V(4)*(Y-Y(4)))/W(4)%求點(diǎn)X(4),Y(4),Z(4)的切平面surf(X,Y,Z); %在曲面上畫法線holdon;surf(X,Y,Z1)；holdon在上面點(diǎn)處畫切平面plot3(X(4),Y(4),Z(4),'r.','MarkerSize',30)%將該點(diǎn)突出放大并用紅色表示U(1:3)=NaN;U(5:end)=NaNholdon;quiver3(X,Y,Z,U,V,W,-4)%作該點(diǎn)處的法線，長度為4,反方向axisequal%該命令是必須的，使得視覺正確5.3.3數(shù)值積分對于形如的積分積分的一般命令為：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表示被積函數(shù)，系數(shù)a，b為積分的上下限，tol為用戶給定的誤差限，如不給出則系統(tǒng)內(nèi)定的值為0.0000001?！纠?.3.6】以下是一個比較著名的例子，計(jì)算積分Q=quad('1./(1+25*x.^2)',-1,1)

ans=0.5494也可以使用內(nèi)聯(lián)函數(shù)定義被積函數(shù)F=inline('1./(1+25*x.^2)');%定義一個內(nèi)聯(lián)函數(shù)Q=quad(F,-1，1);ans=0.5494【例5.3.7】對空間曲線求積分，設(shè)空間曲線的參數(shù)方程為：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)計(jì)算結(jié)果為：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求函數(shù)的根函數(shù)求根或稱求函數(shù)的根，是初等數(shù)學(xué)就學(xué)習(xí)過的。然而對于復(fù)雜的函數(shù)、方程組求根的運(yùn)算是用搜索的方法根據(jù)函數(shù)的特性逐步接近某個根，在使用MATLAB命令時(shí)應(yīng)該對函數(shù)有所了解，給出某個根附近的初始點(diǎn)或根所在的某一個區(qū)間。主要求根命令為：roots：求多項(xiàng)式的根fzero：求一般一元函數(shù)的根fsolve：求非線性方程組的根5.4.1求多項(xiàng)式和一般一元函數(shù)的根【例5.4.1】求多項(xiàng)式的根clear;p=[1,-6,-72,-27]r=roots(p)p=1-6-72-27r=12.1229

-5.7345

-0.3884求一元函數(shù)的0點(diǎn)命令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：被求解的函數(shù)，可用inline或字符串定義x0：為初初始值，可以是一個數(shù)，也可以是一個區(qū)間fval：目標(biāo)函數(shù)值options：輸出顯示項(xiàng)，可選擇display：'off':沒有輸出'iter'：將每一步的迭代結(jié)果都輸出來'notify'：只有當(dāng)?shù)皇諗繒r(shí)才輸出警告信息（系統(tǒng)內(nèi)定值）TolX：迭代終止精度?！纠?.4.2】在[1，2]范圍sinx^2的根并顯示詳細(xì)迭代步驟。options=optimset('Display','iter')x=fzero('sin(x*x)',[12],options)結(jié)果為：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點(diǎn)，可能沒有，可能有多個0點(diǎn)，甚至可能有無窮的0點(diǎn)。因此找不到一個通解，我們可以通過人機(jī)對話，將函數(shù)顯示在屏幕上，在可能的點(diǎn)上用鼠標(biāo)采樣得到初始點(diǎn)，然后讓計(jì)算機(jī)根據(jù)這些點(diǎn)自動分別地求出每一個根。【例5.4.3】編制程序，對任何函數(shù)求根，在運(yùn)行該程序時(shí)將顯示函數(shù)與零點(diǎn)的圖形，用戶在可能的零點(diǎn)附近用鼠標(biāo)點(diǎn)擊采點(diǎn)，程序通過循環(huán)計(jì)算每個采集點(diǎn)的根，并輸出結(jié)果。以函數(shù)

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

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

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

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

%進(jìn)行搜索Optimizationterminatedsuccessfully:thecurrentxsatisfiestheterminationcriteriausingOPTIONS.TolXof1.000000e-004andF(X)satisfiestheconvergencecriteriausingOPTIONS.TolFunof1.000000e-004最優(yōu)點(diǎn)為：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：目標(biāo)函數(shù)系數(shù)向量A：不等式約束矩陣b：不等式約束右邊項(xiàng)系數(shù)向量Aeq：等式約束矩陣beq：等式約束右邊項(xiàng)系數(shù)向量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ī)劃類似，只是其目標(biāo)函數(shù)具有變量的二次項(xiàng)。其模型為：(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為變量二次型矩陣，其他參數(shù)同線性規(guī)劃類似【例5.5.4】求以下二次規(guī)劃的最優(yōu)值與最優(yōu)解。目標(biāo)函數(shù)其中%首先數(shù)據(jù)準(zhǔn)備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)計(jì)算結(jié)果為：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)化問題對于非線性函數(shù)的求最優(yōu)化問題情況較為復(fù)雜，沒有一個固定的方法。非線性函數(shù)的求最優(yōu)化模型為：其中f(x)，c(x)，ceq(x)，為非線形函數(shù)。求非線性函數(shù)的求最優(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為非線性函數(shù)?！纠?.5.5】求以下非線性函數(shù)的最小值及最優(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)計(jì)算結(jié)果為：x=24.0012.0012.0fval=-34565.5.4多目標(biāo)線性規(guī)劃在實(shí)際問題中目標(biāo)函數(shù)可能有多個，例如一個大型的工程其最終的目標(biāo)是多方面的，既要考慮工程的成本盡量少，又要考慮企業(yè)投產(chǎn)帶來的生態(tài)和社會效應(yīng)等多方面的因素。這就是一個多目標(biāo)的規(guī)劃問題。多目標(biāo)規(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(...)其中g(shù)oal：為各目標(biāo)函數(shù)要達(dá)到的上界，goal為向量，維數(shù)是目標(biāo)的個數(shù)。weight：為搜索步長權(quán)重?！纠?.5.6】求以下多目標(biāo)規(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,[])計(jì)算結(jié)果為：最優(yōu)解為：x=1.00000.3333最終多目標(biāo)值：fval=1.6667-1.33335.6傅立葉分析與快速傅立葉變換傅立葉分析在連續(xù)和離散脈沖信號處理上是非常有用的工具，它可以將輸入的波長進(jìn)行變換，分離出各種不同的頻率或者提取輸入波長不同的數(shù)字特征。一般情況下，我們都是對連續(xù)的波進(jìn)行離散化處理，得到一個離散的向量進(jìn)行傅立葉變換，簡記為（DFT）。而快速傅立葉變換（FTT）是處理離散傅立葉變換最強(qiáng)有力的工具，非常適合處理諸如時(shí)間序列、信號和圖象序列數(shù)據(jù)的分析，可以對信號進(jìn)行濾波處理、變換、譜分析和估計(jì)等等?？梢詫⒏盗⑷~分析看成是一個過濾器。例如：

x1,x2,…,xN

X1,X2,…,XN

X1,X2,…,XN

x1,x2,…,Xn傅立葉變換公式為：DFTIDFT傅立葉的逆變換公式為：如果x（n）為實(shí)數(shù)，則可寫成這里【例5.6.1】這是時(shí)間序列課程中一個著名的例子，觀測太陽黑子變化規(guī)律。我們已經(jīng)有300多年有關(guān)太陽黑子的數(shù)據(jù)，稱為著名的Wolfernumber。并且我們也已經(jīng)知道太陽黑子的變化周期是11年，本例用fft方法對該時(shí)間序列進(jìn)行變換，并進(jìn)行周期特征的提取。loadsunspot.dat%讀入太陽黑子數(shù)據(jù)year=sunspot(:,1);wolfer=sunspot(:,2);plot(year,wolfer)title('300年的太陽黑子數(shù)據(jù)圖')我們可以看出太陽黑子是呈周期性變化的，我們初步觀察它的周期是在10~12年之間。為了得出更精確的周期，我們進(jìn)行序列的特征提取，即進(jìn)行fft變換。Y=fft(wolfer)%對序列進(jìn)行快速傅立葉變換通過變換Y的元素已為復(fù)數(shù)，我們來尋找最大周期。首先求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最高點(diǎn)處周期的值大約為11年。為精確起見，我們來搜索最高點(diǎn)處的周期值。[mp,index]=max(power);period(index)ans=11.0769即周期為11.0769年5.7求解常微分方程

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