第5章 數(shù)值分析

1、第5章 數(shù)值分析在數(shù)學(xué)上有很多問題不能得到滿意的結(jié)果。例如對于一個(gè)函數(shù)積分或微分，不能得到解析解時(shí)，可以借助計(jì)算機(jī)的數(shù)值計(jì)算功能，在數(shù)值上求得近似解。我們稱這種方法為數(shù)值分析。插值法是古老而實(shí)用的數(shù)值方法。插值法是函數(shù)逼近的重要方法。在生產(chǎn)和科學(xué)實(shí)踐中，函數(shù)f (x)有時(shí)不能直接寫出表達(dá)式，而只能給出函數(shù)在若干個(gè)點(diǎn)處對應(yīng)的函數(shù)值或?qū)?shù)值。當(dāng)要求知道觀測點(diǎn)之外的函數(shù)值時(shí)，需要估計(jì)函數(shù)值在任意點(diǎn)上的值（之所以叫估計(jì)，是因?yàn)楸徊搴瘮?shù)的未知性）。如何根據(jù)有限的觀測點(diǎn)的值，構(gòu)造出一個(gè)比較簡單的函數(shù)(x)，使函數(shù)在觀測點(diǎn)的值等于對應(yīng)的數(shù)值或?qū)?shù)值。可以用函數(shù)(x)在任意點(diǎn)x處的值來估計(jì)被插函數(shù)f (x)在

2、x點(diǎn)的值。尋找這樣的函數(shù)(x)，其方法很多，如：1. (x)可以是一個(gè)代數(shù)多項(xiàng)式，或是三角多項(xiàng)式也可以是有理分式；2. (x)可以是任意光滑（任意階導(dǎo)數(shù)連續(xù)）的函數(shù)或是分段函數(shù)。函數(shù)類的不同，自然有不同的逼近效果。在許多應(yīng)用中，我們通常要用一個(gè)解析函數(shù)（一、二元函數(shù)）來描述數(shù)據(jù)（通常是實(shí)踐中的測量數(shù)據(jù)）。分兩種情況考慮：1. 測量值準(zhǔn)確，沒有誤差；2. 測量值與真實(shí)值有誤差，這時(shí)有兩種方法： 插值或曲線擬合； 回歸分析。當(dāng)假定數(shù)據(jù)正確時(shí)，一般用插值法，否則用曲線擬合。5.1級數(shù)求和在數(shù)學(xué)中，級數(shù)求和有數(shù)項(xiàng)級數(shù)和函數(shù)項(xiàng)級數(shù)之分，且并不是所有的級數(shù)都是收斂的，但在Matlab中沒有考慮，只要給出數(shù)

3、列的通式，用符號工具箱（Symbolic Math Toolbox）中的和命令symsum對通式操作即可，其常用形式為：S=symsum(s)S=symsum(s,v)S=symsum(s,v,a,b)參數(shù)說明：（1）s：數(shù)列的通項(xiàng)式；（2）v：通式s中的變量，求和時(shí)將對v 從1求至v-1，若不指定v，則對缺省變量求和，如k從0至k-1求和；（3）a,b：對變量從a至b求和，b可以為無窮大，即求和至無窮項(xiàng)；（4）S：求和結(jié)果。對于一些求不出的和數(shù)，Matlab會(huì)給出求和形式。例5-1 求下列數(shù)列的和（1）S=(2) S=(3) S=解：>> syms n a q>> s

4、1=1/(2*n-1)2;>> s2=a*qn;>> s3=1/(n*(n+1)*(n+2)2;>> S1=symsum(s1) S1 = -1/4*Psi(1,n-1/2)>> S2=symsum(s2,n) S2 = a*qn/(q-1)>> S3=symsum(s3,n,1,inf) S3 = -39/16+1/4*pi2不過應(yīng)注意的一點(diǎn)是，并非所有在數(shù)學(xué)上能求和的級數(shù)，在Matlab中都能計(jì)算。5.2 插值擬合用一條曲線去擬合一些離散數(shù)據(jù)，會(huì)面臨兩個(gè)問題：1. 曲線與數(shù)據(jù)之間的擬合好壞標(biāo)準(zhǔn)。2. 用何種函數(shù)類來擬合數(shù)據(jù)。下面分

5、別來回答這兩個(gè)問題。對于第一個(gè)問題，當(dāng)然要求曲線越接近每個(gè)數(shù)據(jù)點(diǎn)越好，但有時(shí)很難做到。一般只要求曲線到數(shù)據(jù)點(diǎn)之間的距離的平方和越小越好，即平常所講的最小二乘曲線擬合。對于第二個(gè)問題，理論上應(yīng)該可用任意的函數(shù)類來進(jìn)行擬合。但有時(shí)需要計(jì)算數(shù)據(jù)的其它值，這時(shí)對擬合函數(shù)要求有很好的解析性，如可導(dǎo)，可積分等。而多項(xiàng)式函數(shù)都具有上述性質(zhì)。5.2.1 一維擬合在Matlab中，一般選用多項(xiàng)式作為擬合函數(shù)。多項(xiàng)式擬合的命令有：1. p = polyfit (x, y, n)參數(shù)說明：·x, y：測量數(shù)據(jù)的橫縱坐標(biāo)向量，n為多項(xiàng)式的次數(shù)（若有N個(gè)數(shù)據(jù)點(diǎn)，則n最大為N-1，因?yàn)镹-1次多項(xiàng)式有N個(gè)參數(shù)，

6、而N個(gè)數(shù)據(jù)點(diǎn)恰好確定N個(gè)等式。注意：并非n越大越好?。?；·p：擬合多項(xiàng)式的系數(shù)向量（按降冪排列）。2. Y = polyval (p, x)參數(shù)說明：·p：為多項(xiàng)式的系數(shù)向量（p可能是一般多項(xiàng)式，或者是擬合多項(xiàng)式）；·x：為需要計(jì)算的自變量向量（x中的值可以落在擬合區(qū)間之外?。?；·Y：為多項(xiàng)式在x處的值。例5-2 在某次工程測量中得到如下數(shù)據(jù)：X1520253035404550Y07442267034163549785502095022650230對以上數(shù)據(jù)用一次數(shù)合適的多項(xiàng)式進(jìn)行擬合，并畫圖比較。解：先把數(shù)據(jù)點(diǎn)描出來：>> x=15:5:5

7、0;>> y=0 7442 26703 41635 49785 50209 50226 50230;>> subplot(2,2,1),plot(x,y,'*') %在plot命令前加上subplot (m, n, p)表示將窗口劃分為m×n個(gè)區(qū)域，而下一個(gè)plot命令則會(huì)繪圖于從左至右的第p個(gè)區(qū)域。再分別用2次，3次，5次多項(xiàng)式對x， y進(jìn)行擬合，分別畫出圖形，與原數(shù)據(jù)點(diǎn)進(jìn)行比較：>>p2=polyfit(x,y,2)p2 = 1.0e+004 * -0.0066 0.5845 -7.6669>> p3=polyfit

8、(x,y,3)p3 = 1.0e+004 * -0.0001 0.0035 0.2800 -4.8880>> p5=polyfit(x,y,5)Warning: Polynomial is badly conditioned. Remove repeated data points or try centering and scaling as described in HELP POLYFIT.> In E:matlab6p1toolboxmatlabpolyfunpolyfit.m at line 74p5 = 1.0e+005 * -0.0000 0.0000 -0.0

9、015 0.0526 -0.8082 4.5106>> xi=15:5:50;>> y2i=polyval(p2,xi);>> y3i=polyval(p3,xi);>> y5i=polyval(p5,xi);>> subplot(2,2,2),plot(x,y,xi,y2i,':')>> subplot(2,2,3),plot(x,y,xi,y3i,'-')>> subplot(2,2,4),plot(x,y,xi,y5i,'-')從上圖中可以看出，5次多項(xiàng)式的

10、擬合效果是最好的。5.2.2 一維插值Matlab中的一維插值函數(shù)命令interp1的常用形式為：t = interp1(X, Y, X0, method)參數(shù)說明：·X：原始數(shù)據(jù)的橫坐標(biāo)向量，必須是單調(diào)增加的向量；·Y：原始數(shù)據(jù)的縱坐標(biāo)向量；·X0：待插值的點(diǎn)的橫坐標(biāo)，可以是標(biāo)量或單調(diào)增加的向量；·t：若X0為標(biāo)量，xi<x0<xi+1，則該方式用的是直線插值。即認(rèn)為被插函數(shù)在點(diǎn)X0處的值落在過點(diǎn) (xi, yi)與(xi+1, yi+1)的直線上，其值為t。若X0為向量，則t也是同維向量；·method：指定插值的算法，取值如下

11、： nearest：最近插值。該方式不進(jìn)行插值，而是找出與點(diǎn)X0最接近的原始數(shù)據(jù)點(diǎn)，再返回其值； linear：線性插值。該方式認(rèn)為相鄰的點(diǎn)之間是線性關(guān)系； cubic：三次插值。該方式認(rèn)為被插函數(shù)在點(diǎn)X0處的值落在過相鄰兩點(diǎn)的三次曲線上； spline：樣條插值。該方式認(rèn)為被插函數(shù)在點(diǎn)X0處的值落在過相鄰兩點(diǎn)的三次樣條曲線上。例5-3 在例5-2中，由于實(shí)際的需要，要估計(jì)x = 22、27、36時(shí)y的值，比較幾種插值方式的差異。（比較圖如下所示）解：>> x=15:5:50;>> y=0 7442 26703 41635 49785 50209 50226 50230

12、;>> X0=22 27 36;>> Y0_n=interp1(x,y,X0,'nearest');>> Y0_c=interp1(x,y,X0,'cubic');>> Y0_s=interp1(x,y,X0,'spline');>> Y0_l=interp1(x,y,X0,'linear');>> p5=polyfit(x,y,5);Warning: Polynomial is badly conditioned. Remove repeated data

13、points or try centering and scaling as described in HELP POLYFIT.> In E:matlab6p1toolboxmatlabpolyfunpolyfit.m at line 74>> xi=15:5:50;>> y5i=polyval(p5,xi);>> plot(x,y,'*',xi,y5i,':',X0,Y0_n,'s',X0,Y0_c,'p',X0,Y0_s,'h',X0,Y0_l,'d'

14、)運(yùn)行結(jié)果為：5.2.3 二維插值在地質(zhì)測量中，通常要用到二維插值。如要繪制一塊矩形土地的地形圖，不可能對該區(qū)域中的每一個(gè)點(diǎn)進(jìn)行測量，只能從中取樣，測量有限個(gè)點(diǎn)的高度，再用插值法確定其它點(diǎn)的高度。在Matlab中二維插值的命令為interp2，常用形式為：Z0 = interp2 (X, Y, Z, X0, Y0, method)參數(shù)說明：·X，Y：原始數(shù)據(jù)的x, y-坐標(biāo)向量，X，Y必須是單調(diào)增加的向量；·Z：原始數(shù)據(jù)的z-坐標(biāo)矩陣，x-y平面上的點(diǎn)(Xi,Yj)與之對應(yīng)的Z的值為Z(i, j)；·X0，Y0：插值點(diǎn)的x, y-坐標(biāo)。若插值點(diǎn)只有一個(gè)，則X0，Y

15、0為標(biāo)量，若插值點(diǎn)有多個(gè)，則為向量；·method：指定插值的算法。取值如下： linear：線性插值。該方式認(rèn)為被插函數(shù)在點(diǎn)(X0, Y0)處的值是在過四個(gè)相鄰點(diǎn)（因?yàn)橥ǔ樾【匦危┑钠嫔希?cubic：三次插值。該方式認(rèn)為被插函數(shù)在點(diǎn)(X0, Y0)處的值落在過某小區(qū)域的三次曲面上； nearest：該方式不進(jìn)行插值。而是找出與點(diǎn)(X0, Y0)最接近的原始數(shù)據(jù)點(diǎn)，再返回其值； spline：樣條插值。該方式認(rèn)為被插函數(shù)在點(diǎn)(X0, Y0)處的值落在過某小區(qū)域的三次樣條曲面上；·Z0：被插函數(shù)在點(diǎn)(X0, Y0)的值。當(dāng)X0與Y0為向量時(shí)，Z0為插值矩陣。例5-4 對

16、山峰函數(shù)peaks進(jìn)行插值，并對不同的插值算法進(jìn)行比較。解：>> x,y,z=peaks(10);>> subplot(2,2,1)>> mesh(x,y,z) %原始數(shù)據(jù)>> xi,yi=meshgrid(-3:0.1:3,-3:.1:3); %生成更細(xì)小的插值點(diǎn)>> zi_l=interp2(x,y,z,xi,yi,'linear'); %用不同的插值算法>> zi_c=interp2(x,y,z,xi,yi,'cubic');>> zi_s=interp2(x,y,z,xi

17、,yi,'spline');>> subplot(2,2,2)>> mesh(xi,yi,zi_l) %線性插值圖>> subplot(2,2,3)>> mesh(xi,yi,zi_c) %三次插值圖>> subplot(2,2,4)>> mesh(xi,yi,zi_s) %樣條插值圖運(yùn)行結(jié)果為：5.3 三次樣條 5.3.1 概述對于給定的離散的測量數(shù)據(jù)（稱為斷點(diǎn)），要尋找一個(gè)三次多項(xiàng)式，y = p (x)，以逼近每對數(shù)據(jù)點(diǎn)間的曲線。過兩點(diǎn)(xi, yi)和(xi+1, yi+1)只能確定一條直線，而通過兩

18、點(diǎn)的三次多項(xiàng)式曲線有無窮多條。為使三次多項(xiàng)式曲線具有唯一性，要增加兩個(gè)條件（因?yàn)槿味囗?xiàng)式有4個(gè)系數(shù)）：1. 三次多項(xiàng)式在點(diǎn)(xi, yi)處有： 2. 三次多項(xiàng)式在點(diǎn)(xi+1, yi+1)處有： 為了使三次多項(xiàng)式具有良好的解析性，還要加上如下條件：3. p (x)在點(diǎn)(xi, yi)處的斜率是連續(xù)的；4. p (x)在點(diǎn)(xi, yi)處的曲率是連續(xù)的；對于第一個(gè)和最后一個(gè)多項(xiàng)式，人為地規(guī)定如下條件：（1）（2）上述兩個(gè)條件稱為非結(jié)點(diǎn)(not-knot )條件?？芍獙?shù)據(jù)擬合的三次樣條函數(shù)p (x)是一個(gè)分段的三次多項(xiàng)式：其中每一段pi (x)都是三次多項(xiàng)式。5.3.2 函數(shù)spline從

19、上面可知，若有N個(gè)數(shù)據(jù)點(diǎn)，要求得一個(gè)三次樣條函數(shù)，則要計(jì)算一個(gè)N*N 的方程組，過程繁瑣。Matlab提供了一個(gè)用于求三次樣條函數(shù)的命令spline，其使用形式為： pp = spline (x, y) y0 = spline (x, y, x0)參數(shù)說明：·x, y：原始的測量數(shù)據(jù)，且x必須是單調(diào)增加的向量；·x0：期望計(jì)算的值；·pp：為包含分段三次多項(xiàng)式的所有信息：有斷點(diǎn)個(gè)數(shù)n，斷點(diǎn)編號，各個(gè)小區(qū)間中三次多項(xiàng)式的系數(shù)；·y0：y0為三次樣條函數(shù)在x0點(diǎn)的值。例5-5 對離散的分布在sin函數(shù)曲線上的數(shù)據(jù)點(diǎn)進(jìn)行樣條插值，畫出插值點(diǎn)圖。解：>&g

20、t; x=0:20;y=sin(x);>> xx=0:.25:20;>> yy=spline(x,y,xx);>> plot(x,y,'o',xx,yy)5.3.3 函數(shù)ppval和unmkpp對于spline的第一種調(diào)用形式，得到的是三次樣條的有關(guān)信息，若要計(jì)算三次樣條在其它點(diǎn)x0的值，或是想對三次樣條的信息有進(jìn)一步的了解，則有命令：y0 = ppval (pp, x0) %相當(dāng)于spline第二種調(diào)用形式breaks, coefs, npolys, ncoefs = unmkpp (pp)參數(shù)說明：·pp：包含樣條函數(shù)所有系數(shù)的

21、多項(xiàng)式；·x0, y0：x0為要計(jì)算的數(shù)值點(diǎn)，y0為相應(yīng)的函數(shù)值；·breaks：斷點(diǎn)數(shù)，即數(shù)據(jù)點(diǎn)數(shù)；·coefs：系數(shù)矩陣。它的第i行是第i個(gè)多項(xiàng)式的系數(shù)；·npolys：多項(xiàng)式的分段數(shù)；·ncoefs：每個(gè)多項(xiàng)式系數(shù)的數(shù)目，一般為4。例5-6 在某次測量中得到如下數(shù)據(jù)：X2.342.983.354.035.006.777.518.63Y23.420.124.625.722.926.829.621.7試用樣條函數(shù)估計(jì)X在其它點(diǎn)3.46, 4.55, 6.03處的值，順便查看樣條函數(shù)的相關(guān)系數(shù)。解：>> X=2.34 2.98 3.

22、35 4.03 5.00 6.77 7.51 8.63;>> Y=23.4 20.1 24.6 25.7 22.9 26.8 29.6 21.7;>> X0=3.46 4.55 6.03;>> pp=spline(X,Y);>> Y0=ppval(pp,X0)Y0 = 25.4466 24.0345 23.8785>> breaks,coefs,npolys,ncoefs=unmkpp(pp)breaks = Columns 1 through 7 2.3400 2.9800 3.3500 4.0300 5.0000 6.7700 7

23、.5100 Column 8 8.6300coefs = -24.3286 57.2891 -31.8563 23.4000 -24.3286 10.5782 11.5788 20.1000 7.2940 -16.4265 9.4149 24.6000 1.5100 -1.5468 -2.8069 25.7000 -0.4120 2.8472 -1.5455 22.9000 -2.4947 0.6597 4.6617 26.8000 -2.4947 -4.8785 1.5398 29.6000npolys = 7ncoefs = 45.4 方程（組）的迭代求解5.4.1 線性方程組的迭代求解設(shè)方程組為 AX = b 對于一般條件下的矩陣A，要用迭代法計(jì)算的解比較繁瑣，我們只在A滿足一定條件下，利用迭代法求的近似解。若A滿足如下條件：1. A為n階方陣，且|A|0；2. A的所有對角線元素aii0則有唯一解。把改寫成迭代形式： X = BX+g其中：B = E-diag (aii-1)A，g = diag (aii-1) b，E為單位矩陣。設(shè)X(0)是一個(gè)任意給定的初始迭代向量，由 X (m)= BX(m-1)+g (m = 1, 2, )得向量序列：X(1)，X(2)，對于