第5講 MATLAB多項(xiàng)式及插值.ppt

1、第5講 MATLAB數(shù)據(jù)分析與多項(xiàng)式計(jì)算,5.1 數(shù)據(jù)統(tǒng)計(jì)處理 5.2 多項(xiàng)式計(jì)算 5.3 插值運(yùn)算 5.4 多項(xiàng)式擬合,5.1 數(shù)據(jù)統(tǒng)計(jì)處理 5.1.1 最大值和最小值,MATLAB提供的求數(shù)據(jù)序列的最大值和最小值的函數(shù)分別為max和min，兩個(gè)函數(shù)的調(diào)用格式和操作過(guò)程類(lèi)似。,1求向量的最大值和最小值 求一個(gè)向量X的最大值的函數(shù)有兩種調(diào)用格式，分別是： (1) y=max(X)：返回向量X的最大值存入y，如果X中包含復(fù)數(shù)元素，則按模取最大值。,(2) y,I=max(X)：返回向量X的最大值存入y，最大值的序號(hào)存入I，如果X中包含復(fù)數(shù)元素，則按模取最大值。 求向量X的最小值的函數(shù)是min(X

2、)，用法和max(X)完全相同。,例5-1 求向量x的最大值。 命令如下： x=-43,72,9,16,23,47; y=max(x) %求向量x中的最大值 y,l=max(x) %求向量x中的最大值及其該元素的位置,2求矩陣的最大值和最小值,求矩陣A的最大值的函數(shù)有3種調(diào)用格式，分別是：,(1) max(A)：返回一個(gè)行向量，向量的第i個(gè)元素是矩陣A的第i列上的最大值。 (2) Y,U=max(A)：返回行向量Y和U，Y向量記錄A的每列的最大值，U向量記錄每列最大值的行號(hào)。 (3) max(A,dim)：dim取1或2。dim取1時(shí)，該函數(shù)和max(A)完全相同；dim取2時(shí)，該函數(shù)返回一個(gè)

3、列向量，其第i個(gè)元素是A矩陣的第i行上的最大值。,求最小值的函數(shù)是min，其用法和max完全相同。,例5-2 分別求34矩陣x中各列和各行元素中的最大值，并求整個(gè)矩陣的最大值和最小值。,3兩個(gè)向量或矩陣對(duì)應(yīng)元素的比較,函數(shù)max和min還能對(duì)兩個(gè)同型的向量或矩陣進(jìn)行比較，調(diào)用格式為：,(1) U=max(A,B)：A,B是兩個(gè)同型的向量或矩陣，結(jié)果U是與A,B同型的向量或矩陣，U的每個(gè)元素等于A,B對(duì)應(yīng)元素的較大者。 (2) U=max(A,n)：n是一個(gè)標(biāo)量，結(jié)果U是與A同型的向量或矩陣，U的每個(gè)元素等于A對(duì)應(yīng)元素和n中的較大者。 min函數(shù)的用法和max完全相同。,例5-3 求兩個(gè)23矩陣

4、x, y所有同一位置上的較大元素構(gòu)成的新矩陣p。,5.1.2 求和與求積,數(shù)據(jù)序列求和與求積的函數(shù)是sum和prod，其使用方法類(lèi)似。設(shè)X是一個(gè)向量，A是一個(gè)矩陣，函數(shù)的調(diào)用格式為：,sum(X)：返回向量X各元素的和。 prod(X)：返回向量X各元素的乘積。 sum(A)：返回一個(gè)行向量，其第i個(gè)元素是A的第i列的元素和。 prod(A)：返回一個(gè)行向量，其第i個(gè)元素是A的第i列的元素乘積。 sum(A,dim)：當(dāng)dim為1時(shí)，該函數(shù)等同于sum(A)；當(dāng)dim為2時(shí)，返回一個(gè)列向量，其第i個(gè)元素是A的第i行的各元素之和。 prod(A,dim)：當(dāng)dim為1時(shí)，該函數(shù)等同于prod(A

5、)；當(dāng)dim為2時(shí)，返回一個(gè)列向量，其第i個(gè)元素是A的第i行的各元素乘積。,例5-4 求矩陣A的每行元素，每列元素的乘積和全部元素的乘積。,5.1.3 平均值和中值,求數(shù)據(jù)序列平均值的函數(shù)是mean，求數(shù)據(jù)序列中值的函數(shù)是median。兩個(gè)函數(shù)的調(diào)用格式為：,mean(X)：返回向量X的算術(shù)平均值。 median(X)：返回向量X的中值。 mean(A)：返回一個(gè)行向量，其第i個(gè)元素是A的第i列的算術(shù)平均值。 median(A)：返回一個(gè)行向量，其第i個(gè)元素是A的第i列的中值。 mean(A,dim)：當(dāng)dim為1時(shí)，該函數(shù)等同于mean(A)；當(dāng)dim為2時(shí)，返回一個(gè)列向量，其第i個(gè)元素是A

6、的第i行的算術(shù)平均值。 median(A,dim)：當(dāng)dim為1時(shí)，該函數(shù)等同于median(A)；當(dāng)dim為2時(shí)，返回一個(gè)列向量，其第i個(gè)元素是A的第i行的中值。,例5-5 別求4元素向量x的平均值和中值。,5.1.4 排序,MATLAB中對(duì)向量X是排序函數(shù)是sort(X)，函數(shù)返回一個(gè)對(duì)X中的元素按升序排列的新向量。,sort函數(shù)也可以對(duì)矩陣A的各列或各行重新排序，其調(diào)用格式為： Y,I=sort(A,dim) 其中dim指明對(duì)A的列還是行進(jìn)行排序。若dim=1，則按列排；若dim=2，則按行排。Y是排序后的矩陣，而I記錄Y中的元素在A中位置。,一個(gè)多項(xiàng)式的冪級(jí)數(shù)形式可表示為： 同時(shí)也可表

7、為嵌套形式 或因子形式 N階多項(xiàng)式n個(gè)根，其中包含重根和復(fù)根。若多項(xiàng)式所有系數(shù)均為實(shí)數(shù)，則全部復(fù)根都將以共軛對(duì)的形式出現(xiàn),5.2 多項(xiàng)式運(yùn)算,1、多項(xiàng)式表示,2、相關(guān)函數(shù),冪系數(shù):在MATLAB里，多項(xiàng)式用行向量表示，其元素為多項(xiàng)式的系數(shù)，并從左至右按降冪排列。 例： 被表示為 p=2 1 4 5 poly2sym(p) ans = 2*x3+x2+4*x+5 Roots: 多項(xiàng)式的零點(diǎn)可用命令roots求的。 例： r=roots(p) 得到 r = 0.2500 + 1.5612i 0.2500 - 1.5612i -1.0000 所有零點(diǎn)由一個(gè)列向量給出。,Poly: 由零點(diǎn)可得原始多項(xiàng)

8、式的各系數(shù)，但可能相差一個(gè)常數(shù)倍。 例：根據(jù)上例： poly(r) ans = 1.0000 0.5000 2.0000 2.5000 注意：若存在重根，這種轉(zhuǎn)換可能會(huì)降低精度。 例： 舍入誤差的影響，與計(jì)算精度有關(guān)。,polyval: 可用命令polyval計(jì)算多項(xiàng)式的值。 例： 計(jì)算y(2.5) c=3,-7,2,1,1； xi=2.5; yi=polyval(c,xi) yi = 23.8125,如果xi是含有多個(gè)橫坐標(biāo)值的數(shù)組，則yi也為與xi長(zhǎng)度相同的向量。 c=3,-7,2,1,1; xi=2.5,3; yi=polyval(c,xi) yi = 23.8125 76.0000,5

9、.3插值運(yùn)算,1、Lagrange插值,方法介紹 對(duì)給定的n個(gè)插值點(diǎn) 及對(duì)應(yīng)的函數(shù)值 ，利用n次Lagrange插值多項(xiàng)式，則對(duì)插值區(qū)間內(nèi)任意x的函數(shù)值y可通過(guò)下式求的： MATLAB中沒(méi)有直接實(shí)現(xiàn)拉格朗日算法的函數(shù)，我們已經(jīng)介紹過(guò)該函數(shù)的書(shū)寫(xiě)：,function y=lagrange (a,b,x) y=0; for i=1:length(a) l=1; for j=1:length(b) if j=i l=l; else l=l.*(x-a(j)/(a(i)-a(j); end end y=y+l*b(i); end,算例：給出f(x)=ln(x)的數(shù)值表，用Lagrange計(jì)算ln(0.

10、54) 的近似值。 x=0.4:0.1:0.8； y=-0.916291,-0.693147,-0.510826,-0.356675,-0.223144; lagrange(x,y,0.54) ans = -0.6161 （精確解-0.616143),2、Runge現(xiàn)象和分段插值,問(wèn)題的提出：根據(jù)區(qū)間a,b上給出的節(jié)點(diǎn)做插值多項(xiàng)式p(x)的近似值，一般總認(rèn)為p(x)的次數(shù)越高則逼近f(x)的精度就越好，但事實(shí)并非如此。 反例： 在區(qū)間-5,5上的各階導(dǎo)數(shù)存在，但在此區(qū)間上取n個(gè)節(jié)點(diǎn)所構(gòu)成的Lagrange插值多項(xiàng)式在全區(qū)間內(nèi)并非都收斂。 取n=10,用Lagrange插值法進(jìn)行插值計(jì)算。, x

11、=-5:1:5; y=1./(1+x.2); x0=-5:0.1:5; y0=lagrange(x,y,x0); y1=1./(1+x0.2); 繪制圖形 plot(x0,y0,-r) 插值曲線(xiàn) hold on plot(x0,y1,-b) 原曲線(xiàn) 為解決Rung問(wèn)題，引入分段插值。,算法分析：所謂分段插值就是通過(guò)插值點(diǎn)用折線(xiàn)或低次曲線(xiàn)連接起來(lái)逼近原曲線(xiàn)。 MATLAB實(shí)現(xiàn) 可調(diào)用內(nèi)部函數(shù)。 命令 interp1,功能 ： 一維數(shù)據(jù)插值（表格查找）。該命令對(duì)數(shù)據(jù)點(diǎn)之間計(jì)算內(nèi)插值。它找出一元函數(shù)f(x)在中間點(diǎn)的數(shù)值。其中函數(shù)f(x)由所給數(shù)據(jù)決定。,格式：yi = interp1(x,Y,xi

12、,method) %用指定的算法計(jì)算插值： nearest：最近鄰點(diǎn)插值，直接完成計(jì)算； linear：線(xiàn)性插值（缺省方式），直接完成計(jì)算； spline：三次樣條函數(shù)插值。 cubic： 分段三次Hermite插值。,t = 1900:10:1990; p = 75.995 91.972 105.711 123.203 131.669. 150.697 179.323 203.212 226.505 249.633; 對(duì)應(yīng)于美國(guó)從1900年到1990年的每10年的人口數(shù)，求1975年的人口。由此推斷美國(guó)1900年到2000年每一年的人口數(shù)，并畫(huà)出圖形。,與1階拉格朗日算法比較：,lagran

13、d(1970 1980,203.212 226.505,1975) ans = 214.8585,推廣到多個(gè)點(diǎn)計(jì)算：, t = 1900:10:1990; p = 75.995 91.972 105.711 123.203 131.669. 150.697 179.323 203.212 226.505 249.633; x=1965 1975 x = 1965 1975 y=interp1(t,p,x) y = 191.2675 214.8585,lagrand(1960 1970,179.323 203.212,1965) ans = 191.2675 lagrand(1970 1980,

14、203.212 226.505,1975) ans = 214.8585,例, x0=-1+2*0:10/10; y0=1./(1+25*x0.2); x=-1:.01:1; y=lagrange(x0,y0,x); % Lagrange 插值 ya=1./(1+25*x.2); plot(x,ya,x,y,:), y1=interp1(x0,y0,x); plot(x,ya,x,y1);,用插值的方法對(duì)一函數(shù)進(jìn)行近似,要求所得到的插值多項(xiàng)式經(jīng)過(guò)已知插值節(jié)點(diǎn);在n比較大的情況下,插值多項(xiàng)式往往是高次多項(xiàng)式,這也就容易出現(xiàn)振蕩現(xiàn)象（龍格現(xiàn)象），即雖然在插值節(jié)點(diǎn)上沒(méi)有誤差,但在插值節(jié)點(diǎn)之外插值誤差

15、變得很大,從“整體”上看,插值逼近效果將變得“很差”。,所謂數(shù)據(jù)擬合是求一個(gè)簡(jiǎn)單的函數(shù),例如是一個(gè)低次多項(xiàng)式,不要求通過(guò)已知的這些點(diǎn),而是要求在整體上“盡量好”的逼近原函數(shù)。這時(shí),在每個(gè)已知點(diǎn)上就會(huì)有誤差,數(shù)據(jù)擬合就是從整體上使誤差,盡量的小一些。,5.4 數(shù)據(jù)最小二乘擬合,多項(xiàng)式最小二乘擬合,n次多項(xiàng)式：,下面我們練習(xí)一個(gè)具體例題,某種鋁合金的含鋁量為x%，其溶解溫度為y攝氏度，由試驗(yàn)測(cè)得的x與y的數(shù)據(jù)表如下，試用最小二乘算法建立x,y的經(jīng)驗(yàn)公式。,（1）根據(jù)所給出數(shù)據(jù)，畫(huà)出圖形，觀察數(shù)據(jù)關(guān)系,（2）通過(guò)圖形，我們可以確定該曲線(xiàn)可以通過(guò)一次方程y=ax+b的形式來(lái)進(jìn)行擬合?？梢源_定，n=1，

16、L=6。,（3）建立法方程組。,帶入具體數(shù)據(jù), A=6 396.9; 396.9 28365.28; B=1458;101176.3; X=inv(A)*B X = 94.7479 2.2412,則 a=94.7479，b=2.2414 擬合方程為 y=94.7479 + 2.2414x,畫(huà)出比較圖形：,x=36.9 46.7 63.7 77.8 84 87.5; y=181 197 235 270 283 292; plot(x,y,x,y,o) hold on x1=35:5:90; y1=2.2412*x1+ 94.7479; plot(x1,y1,r),而該擬合公式在matlab中有內(nèi)

17、部函數(shù)：多項(xiàng)式擬合MATLAB命令：polyfit格式：p=polyfit(x,y,n),其中，xy為原始樣本點(diǎn)構(gòu)成的向量 n為選定的多項(xiàng)式階數(shù) P為多項(xiàng)式系數(shù)（降冪排列）,則上例采用命令求解：, x=36.9 46.7 63.7 77.8 84 87.5; y=181 197 235 270 283 292; a=polyfit(x,y,1) a = 2.2337 95.3524, x0=0:.1:1; y0=(x0.2-3*x0+5).*exp(-5*x0).*sin(x0); p3=polyfit(x0,y0,3); %先進(jìn)行三次擬合 p3 = 2.8400 -4.7898 1.9432 0.0598 多項(xiàng)式如下： 2.8400*x3-4.7898*x2+1.9432*x+0.0598,已知數(shù)據(jù)點(diǎn)來(lái)自 用多項(xiàng)式擬合方法在不同階次下進(jìn)行擬合。,擬合該數(shù)據(jù)的程序,繪制擬合