第7章計算方法的MATLAB實現(xiàn)講稿

1、第7章 計算方法的MATLAB實現(xiàn)第7章 計算方法的MATLAB實現(xiàn) 計算方法主要研究數(shù)學(xué)問題的數(shù)值解，涉及的內(nèi)容廣。利用MATLAB提供的部分函數(shù)解決，可以實現(xiàn)某些情況下的數(shù)值求解。本章作簡要介紹。7.1 一元非線性方程求解 求解一元非線性方程的方法主要有二分法、割線法、牛頓法等。本節(jié)主要介紹兩個可以求解一元非線性方程的函數(shù)，即fzero函數(shù)和roots函數(shù)。7.1.1 fzero函數(shù)用fzero函數(shù)求一元非線性方程的零點。該函數(shù)的簡單調(diào)用格式為：x=fzero(fun,x0)：如果x0為標(biāo)量，則試圖尋找函數(shù)fun在x0附近的零點。fun參數(shù)是一個M文件函數(shù)或匿名函數(shù)的函數(shù)句柄。函數(shù)fzer

2、o返回的值x靠近fun函數(shù)改變符號的位置。如果搜索失敗，則返回NaN。搜索區(qū)間擴展到發(fā)現(xiàn)inf、NaN或復(fù)數(shù)值時終止搜索。如果x0是一個長度為2的向量，則fzero函數(shù)假設(shè)x0是一個區(qū)間，其中fun(x0(1)與fun(x0(2)異號。如果符號相同，則會出錯。給出相應(yīng)函數(shù)值異號的區(qū)間可以保證fzero函數(shù)返回一個fun函數(shù)改變符號的位置附近的值。x,fval=fzero()：返回解x和解x處目標(biāo)函數(shù)fun的值。x,fval,exitflag=fzero()：還返回一個exitflag值，描述fzero函數(shù)的退出條件。返回值及其描述如下所示：>0，表示函數(shù)找到了零值點x；0，表示沒有發(fā)現(xiàn)零

3、值點。例71求方程的解。程序代碼：clear;clc;f=(x)(x+20)*(x+1)*(x-5)*(x-30);%創(chuàng)建匿名函數(shù)x1=fzero(f,-50)%解一元非線性方程x2=fzero(f,1)x3=fzero(f,6)x4=fzero(f,50)運行結(jié)果：x1 = -20x2 = -1x3 = 5x4 = 307.1.2 roots函數(shù)用roots函數(shù)計算多項式的根，該函數(shù)的調(diào)用格式為：r=roots(c)：返回一個列向量，其元素為多項式c的解。例72求多項式方程的解。程序代碼：clear;clc;P=1 -5 1 -5;%表示多項式 r=roots(P)%求多項式方程的解P1=p

4、oly(r)%由多項式的根返回多項式的系數(shù)運行結(jié)果：r = 5.0000 -0.0000 + 1.0000i -0.0000 - 1.0000iP1 = 1.0000 -5.0000 1.0000 -5.0000注：poly和roots互為逆函數(shù)，函數(shù)poly返回多項式的系數(shù)。7.2 非線性方程組的數(shù)值解法 用MATLAB求解線性方程組的解，在6.3中已介紹了基于矩陣變換的直接解法，除此方法之外，還有很多方法，比如，Jocabi迭代法、Gauss-Seidel迭代法和SOR(超松馳)迭代法等，這里不作詳細(xì)介紹。可參閱有關(guān)文獻。實際工程中得到的數(shù)學(xué)模型往往具有非線性的特點，得到其解析解比較困難。

5、一般采用迭代法求解非線性方程組。比較常見的迭代方法有不動點迭代法、Newton迭代法和擬Newton迭代法等幾種。本節(jié)僅介紹不動點迭代法。設(shè)含有個未知數(shù)和個方程的非線性方程組記為：其中，為由個未知數(shù)構(gòu)成的向量，為由個函數(shù)構(gòu)成的向量。將方程組改寫為便于迭代的等價形式：，由此得出不動點迭代法的迭代公式：。如果得到的序列滿足，則就是的不動點。這樣就可以求出線性方程組的解（近似解）。據(jù)此，編寫不動點迭代法的M文件(文件名為example7_2)：function s=example7_2(x,eps)%用不動點迭代法求非線性方程組的解%x為迭代初值，eps為允許誤差值。if nargin=1 eps=

6、1e-6;elseif nargin<1 error returnendx1=example7_2a(x);% example7_2a是函數(shù)的文件名。while norm(x1-x)>=eps x=x1; x1=example7_2a(x);%循環(huán)迭代ends=x1;return例73用不動點迭代法求方程組的解。解：變形為首先創(chuàng)建函數(shù)代碼（文件名為example7_2a.m）：function y=example7_2a(x)y(1)=(x(1)*x(1)+x(2)*x(2)+9)/(-11);y(2)=(x(1)*x(2)*x(2)+x(1)+8)/6;再作圖形，程序代碼：cle

7、ar;clc;ezplot('x12+11*x1+x22+9',-12,6);hold onezplot('x1*x22+x1-6*x2+8',-12,6);hold offtitle('x12+11*x1+x22+9和x1*x22+x1-6*x2+8的圖形')%修改標(biāo)題運行結(jié)果見圖71：圖71最后調(diào)用函數(shù)example7_2求方程組解：x=example7_2 (0,0)運行結(jié)果：x= -1.0000 1.0000 這樣就可以求得方程組的一個解，從圖71可以看出，方程組還有一個解，為了求得它，還需要對原方程組作另外變形，這里不再作進一步的討論。

8、7.3 插值 插值計算在數(shù)據(jù)擬合和數(shù)據(jù)平滑等方面應(yīng)用普遍。MATLAB提供了用最近鄰插值、線性插值、三次樣條插值、三次插值和FFT插值法進行一維、二維、三維和高維插值的函數(shù)。7.3.1 一維插值MATLAB中有兩種一維插值，即多項式插值和基于FFT的插值。 多項式插值函數(shù)interp1進行一維插值。一維插值是進行數(shù)據(jù)分析和曲線擬合的重要手段。interp1函數(shù)使用多項式技術(shù)，用多項式函數(shù)擬合所提供的數(shù)據(jù)，并計算目標(biāo)插值點上的插值函數(shù)值，其最常用的語法形式是：yi=interp1(x,y,xi,method)x和y為給定數(shù)據(jù)的向量，長度相同。xi為包含要插值的點的構(gòu)成的向量，method是一個可

9、選的字符串，指定一種插值方法，包括：最近鄰插值（method='nearest'）：該方法將插值點的值設(shè)置為已知數(shù)據(jù)點中距離最近的點的值；線性插值（method='linear'）：該方法用線性函數(shù)擬合每對數(shù)據(jù)點，并返回xi處的相關(guān)函數(shù)值；三次樣條插值（method='spline'）：該方法用三次樣條函數(shù)擬合每對數(shù)據(jù)點，用spline函數(shù)在插值處進行三次樣條插值；三次插值（method='pchip'或'cubic'）：該方法用pchip函數(shù)對向量x和y進行分段三次Hermite插值。這幾種方法在速度、內(nèi)存和平滑

10、性等方面有所差別，使用時可以根據(jù)需要進行選擇，包括：最近鄰插值是最快的方法，但是，利用它得到的結(jié)果平滑性最差；線性插值比最近鄰插值要占用更多的內(nèi)存，運行時間略長。與最近鄰法不同，它生成的結(jié)果是連續(xù)的，但是在頂點處有坡度變化；三次樣條插值的運行時間相對來說最長，內(nèi)存消耗比三次插值略少。它生成的結(jié)果平滑性最好。但是，如果輸入數(shù)據(jù)不很均勻，可能會得到意想不到結(jié)果；三次插值需要更多內(nèi)存，而且運行時間比最近鄰插值和線性插值的長。但是，使用此法時，插值數(shù)據(jù)及其導(dǎo)數(shù)都是連續(xù)的。例74 程序代碼：x=1:10;y=1800 778 518 5000 980 588 2799 528 6700 598;x1=1

11、:0.1:10;y1=interp1(x,y,x1,'nearest');%一維最近鄰插值plot(x1,y1)y2=interp1(x,y,x1,'linear');%一維線性插值hold onplot(x1,y2,'color','r')hold off運行結(jié)果見圖72。圖72一維最近鄰插值和線性插值例75 程序代碼：x=1:10;y=1800 778 518 5000 980 588 2799 528 6700 598;plot(x,y);x1=1:0.1:10;hold ony1=interp1(x,y,x1,'s

12、pline');%三次樣條插值plot(x1,y1,'color','r')hold off運行結(jié)果見圖73。圖73三次樣條插值 基于FFT的插值函數(shù)interpft用基于FFT的方法進行一維插值。本方法計算包含周期函數(shù)值的向量的傅里葉變換。然后，它用更多的點計算逆傅里葉變換。該函數(shù)的調(diào)用形式為：y=interpft(x,n)其中，x是一個包含周期函數(shù)值的向量，這些值在等間隔的點上采集。n是樣本大小。7.3.2 二維插值二維插值在圖像處理和可視化方面有著很重要的應(yīng)用。MATLAB用函數(shù)interp2進行二維插值。該函數(shù)的一般形式為：ZI=interp2(

13、X,Y,Z,XI,YI,method)其中，Z是一個矩陣數(shù)組，包含二維函數(shù)的值，X和Y為大小相同的數(shù)組，包含相對于Z的給定值。XI和YI為包含插值點數(shù)據(jù)的矩陣，method表示插值方法，為可選參數(shù)。MATLAB提供了三種不同的插值方法進行二維插值：最近鄰插值（method='nearest'）：該方法用分區(qū)域常數(shù)曲面擬合數(shù)據(jù)，插值點的值是最近點的值；雙線性插值（method='linear'）：該方法用雙線性曲面擬合數(shù)據(jù)點，插值點的值是四個最近的值的組合。本方法是分區(qū)域雙線性的，比雙三次插值法快，并且內(nèi)存消耗更少；雙三次插值（method='cubic&

14、#39;）：該方法用雙三次曲面擬合數(shù)據(jù)點，插值點的值是16個最近點的值的組合。本方法是分區(qū)域三次的，結(jié)果的平滑性比前面兩種的都好。注：所有這些方法都要求X和Y數(shù)據(jù)是單調(diào)的，即從點到點，要么總是遞增的，要么總是遞減的。應(yīng)該用meshgrid函數(shù)準(zhǔn)備這些矩陣。例76 程序代碼：%低分辨率的peaks函數(shù)圖形x,y=meshgrid(-4:1:4);z=peaks(x,y);surf(x,y,z)運行結(jié)果見圖74。圖74低分辨率的peaks函數(shù)圖形例77 程序代碼：x,y=meshgrid(-4:1:4);z=peaks(x,y);xI,yI=meshgrid(-4:0.25:4);zI=inter

15、p2(x,y,z,xI,yI,'nearest');%二維最近鄰插值surf(xI,yI,zI)運行結(jié)果見圖75。圖75二維最近鄰插值例78 程序代碼：x,y=meshgrid(-4:1:4);z=peaks(x,y);xI,yI=meshgrid(-4:0.25:4);zI=interp2(x,y,z,xI,yI,'linear');%二維雙線性插值,'linear'可省略.surf(xI,yI,zI)運行結(jié)果見圖76。圖76二維雙線性插值例79 程序代碼：x,y=meshgrid(-4:1:4);z=peaks(x,y);xI,yI=mesh

16、grid(-4:0.25:4);zI=interp2(x,y,z,xI,yI,'cubic');%二維雙三次插值surf(xI,yI,zI)運行結(jié)果見圖77。圖77二維雙三次插值7.3.3 多維插值MATLAB提供了幾種多維數(shù)據(jù)的插值函數(shù)，如表71中所示。表71多維數(shù)據(jù)的插值函數(shù)函數(shù)描述interp3三維數(shù)據(jù)插值interpn多維數(shù)據(jù)插值ndgrid多維數(shù)據(jù)網(wǎng)格化函數(shù)interp3進行三維插值。計算三維樣本集V中數(shù)據(jù)點之間的值。該函數(shù)的一般形式為：VI=interp3(X,Y,Z,V,XI,YI,ZI,method)其中，X，Y和Z指定數(shù)據(jù)點；V包含與X，Y和Z對應(yīng)的值；XI，

17、YI和ZI為interp3函數(shù)對V中數(shù)據(jù)進行插值的點。對于超出范圍的值，interp3函數(shù)返回NaN。method表示插值方法，為可選參數(shù)。對于三維數(shù)據(jù)，有三種不同的插值方法。最近鄰插值（method='nearest'）：該方法選擇最近點的值；線性插值（method='linear'）：該方法基于最近的8個點進行分區(qū)域線性插值；三次插值（method='cubic'）：該方法基于最近的64個點進行分區(qū)域三次插值。用interpn函數(shù)進行更高維數(shù)據(jù)的插值，該函數(shù)的常用形式為：VI=interpn(X1,X2,X3,V,Y1,Y2,Y3,metho

18、d)其中，Y1,Y2,Y3, 為interpn函數(shù)對V中數(shù)據(jù)進行插值的點。對于超出范圍的值，interpn函數(shù)返回NaN。method表示插值方法，為可選參數(shù)。高維數(shù)據(jù)插值，同樣有最近鄰插值、線性插值和三次插值三種方法。用ndgrid函數(shù)為高維函數(shù)評價和插值生成數(shù)據(jù)數(shù)組。該函數(shù)將一系列輸入向量指定的圖域轉(zhuǎn)換為一系列輸出數(shù)組。第i維是輸入向量xi的拷貝。ndgrid函數(shù)的語法格式為：X1,X2,X3,= ndgrid(x1,x2,x3,)7.4 曲線擬合所謂曲線擬合，就是利用兩個或多個變量的離散點，用平滑曲線來擬合它們之間的關(guān)系。根據(jù)擬合方法的不同，有參數(shù)擬合和非參數(shù)擬合之分。參數(shù)擬合，曲線不要

19、求通過所有點，采用最小二乘法；非參數(shù)擬合，要求曲線通過所有點，采用插值法。由于曲線擬合是數(shù)據(jù)分析最常見的任務(wù)之一，MATLAB提供了多種函數(shù)和工具來進行曲線擬合，另外還有曲線擬合工具箱。7.4.1 最小二乘法最小二乘法通過最小化殘差的平方和來獲得待定系數(shù)的估計。第個數(shù)據(jù)點的殘差定義為測量響應(yīng)值和擬合響應(yīng)值之間的差值，即，殘差的平方和。常見的最小二乘法包括線性最小二乘、加權(quán)線性最小二乘、穩(wěn)健最小二乘和非線性最小二乘等。求解非線性最小二乘問題的Gauss-Newton法和Levenberg-Marquart法是老牌算法。7.4.2 多項式曲線擬合用polyfit函數(shù)計算擬合數(shù)據(jù)集的多項式在最小二乘

20、意義上的系數(shù)，調(diào)用形式為：p=polyfit(x,y,n)x和y是包含要擬合的和數(shù)據(jù)的向量，n是多項式的階次。例710 程序代碼：x=1:10;y=1 3 31 133 381 871 1723 3081 5113 8011;P=polyfit(x,y,4)%多項式曲線擬合，返回多項式的系數(shù)。運行結(jié)果：P = 1.0000 -2.0000 -0.0000 1.0000 1.00007.4.3 相關(guān)工具MATLAB支持用基本擬合界面進行擬合。該擬合界面具有擬合快速，操作簡便的優(yōu)勢。它具有如下功能：使用樣條插值、分段三次艾爾米特插值（PCHIP）或者是1到10階的多項式插值進行數(shù)據(jù)擬合；利用一組數(shù)

21、據(jù)可以同時作多條擬合曲線；可以繪制殘差圖；可以檢查擬合結(jié)果；可以對擬合進行內(nèi)插或外推；用擬合結(jié)果和標(biāo)準(zhǔn)殘差在圖中進行注釋；可以將擬合和計算結(jié)果保存到MATLAB工作空間。下面結(jié)合一個具體的例子加以說明。例711某商店某種產(chǎn)品的銷售量如表72所示。表72某產(chǎn)品銷售量資料年份200120022003200420052006200720082009銷售量（萬件）10.018.025.030.535.038.040.039.538.0請用多項式曲線擬合上述數(shù)據(jù)。按照下面的步驟進行操作。用上述數(shù)據(jù)繪散點圖和線形圖的組合圖（見圖78）；程序代碼：%說明多項式曲線擬合界面的使用方法x=2001:2009;y

22、=10 18 25 30.5 35 38 40 39.5 38;scatter(x,y,'filled','r')%繪散點圖hold onplot(x,y)hold off圖78散點圖和線形圖的組合圖在圖形窗口的“Tools”菜單中單擊“Basic Fitting”菜單選項；兩次單擊“”按鈕。打開的曲線擬合界面“Basic Fitting”對話框如圖79所示。基本擬合界面中各選項的功能包括：Select data下拉式列表框：該列表框中顯示了圖形窗口中圖形用到的所有數(shù)據(jù)集的名稱，在其中選擇要擬合的數(shù)據(jù)。一次只能選擇一組數(shù)據(jù)，但在一組數(shù)據(jù)里可以同時擬合多條曲線。C

23、enter and scale x data單選鈕：選擇此項以后，數(shù)據(jù)中心化為具有0均值，比例化為具有單位標(biāo)準(zhǔn)差。對數(shù)據(jù)進行中心化和比例化，可以提高后面數(shù)值計算的精度。Plot fits：使用本面板，可以用圖形查看當(dāng)前數(shù)據(jù)集的一種或多種擬合結(jié)果。Check to display fits on figure：選擇當(dāng)前數(shù)據(jù)集的擬合類型。有兩種擬合類型可供選擇，即插值和多項式。進行三次樣條插值使用Spline函數(shù)，保形（shape-preserving）插值使用pchip函數(shù)（三次插值）。多項式擬合使用polyfit函數(shù)?？梢赃x擇多種擬合類型。Show equations單選鈕：選此項，在擬合圖形

24、上顯示方程。Plot residuals單選鈕：選此項，顯示擬合曲線的殘差，可用條形圖、散點圖或線形圖顯示。Show norm of residuals單選鈕：選此項，顯示殘差的范數(shù)。殘差的范數(shù)是表示擬合優(yōu)度的一個統(tǒng)計量，值越小，表示擬合程度越高。用norm函數(shù)進行計算，即norm(V,2)，其中V為殘差向量。圖79基本擬合界面Numerical results方框：使用該面板，可以在不繪擬合圖的情況下探察對當(dāng)前數(shù)據(jù)集進行單次擬合的數(shù)值結(jié)果。Fit：選擇擬合當(dāng)前數(shù)據(jù)集的方程。擬合結(jié)果顯示在菜單下面的列表框中。注意，在該菜單中選擇方程并不影響“Plot fits”面板的狀態(tài)。所以。如果試圖在數(shù)據(jù)

25、散點圖中顯示擬合曲線，可能需要在“Plot fits”面板中選擇有關(guān)的核選框。Coefficients and norm of residuals：顯示“Fit”中選擇的方程的計算結(jié)果。Save to workspace：打開一個對話框，使用它將擬合結(jié)果保存到工作空間變量中。Find Y=f(X)：對當(dāng)前的擬合進行內(nèi)插或外推。Enter value(s)：輸入一個MATLAB表達式，進行擬合計算。單擊“Evaluate”按鈕以后計算表達式，結(jié)果顯示在有關(guān)的表格中。當(dāng)前擬合顯示在“Fit”菜單中。Save to workspace：打開一個對話框，使用它將擬合結(jié)果保存到工作空間變量中。Plot

26、evaluated results：選擇此項，計算結(jié)果顯示到圖中的數(shù)據(jù)點上。在基本擬合界面中作以下設(shè)置（見圖710）：用二次多項式擬合數(shù)據(jù)；在擬合圖上顯示方程；將擬合殘差作為條形圖顯示，并將殘差的圖形作為子圖顯示；顯示殘差的范數(shù)。當(dāng)前數(shù)據(jù)集用二次多項式擬合數(shù)據(jù)顯示方程繪殘差圖顯示殘差的范數(shù) 圖710進行選項設(shè)置利用“Plot fits”面板可以可視地探察當(dāng)前數(shù)據(jù)集的多個擬合圖形。為了進行比較，通過選擇合適的核選框來擬合數(shù)據(jù)的其他方程。如果某個方程生成的結(jié)果在數(shù)值上不精確，MATLAB會顯示一個警告信息。此時，應(yīng)該選擇“Center and scale x data” 核選框來改進數(shù)值精度。最后

27、將擬合結(jié)果和殘差顯示在圖711中。圖例顯示了數(shù)據(jù)集和方程的名稱。如果圖例覆蓋了圖形的一部分，可以通過單擊和拖拉操作將它移到其他地方。圖711擬合結(jié)果和殘差圖可以指定一個包含x值的向量，計算這些位置上的擬合值。將該向量輸入到“Evaluate”按鈕左邊的文本框中，然后單擊“Evaluate”按鈕。例如，輸入向量2010:2015后，再單擊“Evaluate”按鈕，即可得到2010年到2015年的銷售量的估計。如圖712所示。選擇“Plot evaluated results”核選框，顯示當(dāng)前圖形中對應(yīng)當(dāng)前數(shù)據(jù)的計算值，如圖713所示圖712顯示x值和對應(yīng)的擬合值圖713擬合結(jié)果和對應(yīng)的殘差圖7.

28、5 數(shù)值微分利用MATLAB函數(shù)，可以實現(xiàn)數(shù)值微分、梯度運算。7.5.1 數(shù)值微分運算5.2一節(jié)中用diff函數(shù)實現(xiàn)了符號微分運算。實際上，diff函數(shù)還可以實現(xiàn)數(shù)值微分運算，即求相鄰元素之間的差值。diff函數(shù)的調(diào)用格式如下：Y=diff(X,n,dim)n和dim均為大于或等于1的整數(shù)， n和dim的默認(rèn)值均為1。當(dāng)dim大于X的維數(shù)時，將返回空數(shù)組；當(dāng)n大于或等于X的第dim維的長度時，也返回空數(shù)組。若X是向量，則可省略dim，diff(X)返回X的一階差分；diff(X,n)(n大于1且小于X的長度)返回X的n階差分(等價于遞歸n次一階差分)。若X是p行q列矩陣，則diff(X)返回X的每列元素的一階差分；diff(X,n)(n大于1且小于p)返回X的每列元素的n階差分；diff(X,n,2)(n大于1且小于q)返回X的每行元素的n階差分。若X是數(shù)組，則diff(X,n,dim)()返回X的第d