1、第三講 數(shù)值運算基礎(chǔ), matlab 具有出色的數(shù)值計算能力，占據(jù)世界上數(shù)值計算軟件的主導(dǎo)地位,數(shù)值運算的功能,多項式運算 線性方程組 數(shù)值統(tǒng)計 線性插值和擬合,matlab語言把多項式表達成一個行向量， 該向量中的元素是按多項式降冪排列的。 f(x)=anxn+an-1xn-1+a0 可用行向量 p=an an-1 a1 a0表示 poly 產(chǎn)生特征多項式系數(shù)向量 特征多項式一定是n+1維的 特征多項式第一個元素一定是1,一、 多項式運算,例: a=1 2 3;4 5 6;7 8 0; p=poly(a) p =1.00 -6.00 -72.00 -27.00 p是多項式p(x)=x3-6x

2、2-72x-27的matlab描述方法 p1=poly2sym(p,x) 顯示數(shù)學(xué)多項式的形式 p1 =x3 6*x2 72*x 27,2. 多項式求值,有兩種求多項式值的函數(shù)：polyval與polyvalm，前者是代數(shù)多項式求值，后者是矩陣多項式求值. (1) polyval函數(shù): y=polyval(p,x) 若x為一數(shù)值，則求多項式在該點的值；若x為向量或矩陣，則對向量或矩陣中的每個元素求其多項式的值。,例3-1 已知多項式x4+8x3-10，分別取x=1.2和一個22矩陣為自變量計算它的值.,eg3-1-2,2. 多項式求值,(2) polyvam函數(shù): y=polyvam(p,x)

3、 本函數(shù)要求x為方陣，它以方陣為自變量求多項式的值. 設(shè)A為方陣，P代表多項式x3-5x2+8，那么polyvalm(P,A)的含義是：A*A*A-5*A*A+ 8*eye(size(A) 而polyval(P,A)的含義是：A.*A.*A-5*A.*A+ 8*ones(size(A),例3-2 仍以多項式x4+8x3-10為例，取一個22矩陣為自變量分別用polyval和polyvalm計算該多項式的值.,eg3-1-2,3.roots 求多項式的根,a=1 2 3;4 5 6;7 8 0;p=poly(a) p = 1.00 -6.00 -72.00 -27.00 %a的特征多項式系數(shù) r

4、=roots(p) r = 12.12 -5.73 -0.39,顯然 r是矩陣a的特征值,當(dāng)然我們可用poly命令將根返回為多項式形式 p2=poly(r) p2 = 1.00 -6.00 -72.00 -27.00 matlab規(guī)定多項式系數(shù)向量用行向量表示，一組根用列向量表示。,根 多項式,例3-3 用求特征值的方法解方程 3x5-7x4+5x2+2x-18=0,eg3-3,命令如下： p=3 -7 0 5 2 -18;A=compan(p) %p的伴隨矩陣x1=eig(A) %求A的特征值 x2=roots(p) %直接求多項式p的零點 即 方程p(x)=0的解,4.conv多項式乘運算

5、,例: a(x)=x2+2x+3; b(x)=4x2+5x+6; c = (x2+2x+3)(4x2+5x+6) a=1 2 3;b=4 5 6; c=conv(a,b)=conv(1 2 3,4 5 6) c = 4.00 13.00 28.00 27.00 18.00 p=poly2sym(c) p = 4*x4 + 13*x3 + 28*x2 + 27*x + 18,5.deconv多項式除運算,a=1 2 3; c = 4.00 13.00 28.00 27.00 18.00 d=deconv(c,a) d =4.00 5.00 6.00,6.多項式微積分,命令格式： polyder(

6、p): 求p的微分 polyder(a,b): 求多項式a,b乘積的微分 p,q=polyder(a,b)：求a/b的微分,分子存入p,分母存入q polyint(p)： 求多項式p的積分 例：a=1 2 3 4 5; poly2sym(a) ans = x4 + 2*x3 + 3*x2 + 4*x + 5 b=polyder(a); poly2sym(b) ans =4*x3 + 6*x2 + 6*x + 4 c=polyint(b); poly2sym(c) ans=x4+2*x3+3*x2+4*x,二、代數(shù)方程組求解,對于方程組ax=b，a 為amn矩陣，有三種情況： 當(dāng)m=n時，此方程

7、稱為“恰定”方程組 當(dāng)mn時，此方程稱為“超定”方程組 當(dāng)mn時，此方程稱為“欠定”方程組 matlab定義的兩種除運算（左除和右除）可以很方便地解上述三種方程,1.恰定方程組的解,方程ax=b(a為非奇異) x=a-1 b 矩陣逆 兩種解: x=inv(a)b 采用求逆運算解方程 x=ab 采用左除運算解方程,方程ax=b a=1 2;2 3;b=8;13; x=inv(a)*b x=ab x = x = 2.00 2.00 3.00 3.00,例: x1+2x2=8 2x1+3x2=13,2.超定方程組的解,方程 ax=b, mn時此時不存在唯一解 方程解 (a a)x=a b x=(a

8、a)-1 a b 求逆法 x=ab matlab用最小二乘法找一 個準(zhǔn)確地基本解。,例: x1+2x2=1 2x1+3x2=2 3x1+4x2=3 a=1 2;2 3;3 4;b=1;2;3; 解1 x=ab 解2 x=inv(aa) a b x = x = 1.00 1.00 0 0.00,3.欠定方程組的解,當(dāng)方程數(shù)少于未知量個數(shù)時,即不定情況,有無窮多個解存在 matlab可求出兩個解： 用除法求的解x是具有最多零元素(即最少元素)的解； 用求逆法的解是具有最小長度或范數(shù)的解，這個解是基于偽逆pinv求得的。,x1+2x2+3x3=1 2x1+3x2+4x3=2 a=1 2 3;2 3

9、4;b=1;2; x=ab x=pinv(a)b x = x = 1.00 0.83 0 0.33 0 -0.17,三、數(shù)據(jù)分析,max 各列最大值 mean 各列平均值 sum 各列求和 std 各列標(biāo)準(zhǔn)差 var 各列方差 sort 各列遞增排序,y=max(X)：返回向量X的最大值存入y，如果X中包含復(fù)數(shù)元素，則按模取最大值。,1. 最大值和最小值(max和min),eg3-8,對矩陣A來說： max(A) 返回一個行向量，向量的第i個元素是矩陣A的第i列上的最大值； 求最小值的函數(shù)min有類似的用法,數(shù)據(jù)序列求和與求積的函數(shù)是sum和prod： sum(A)：返回一個行向量,第i個元素

10、是A的第i列元素和 prod(A)：返回一個行向量,第i個元素是A的第i列的元素乘積,2. 求和與求積,3. 平均值和中值,求數(shù)據(jù)序列平均值與中值的的函數(shù)是mean和median： mean(A)：返回一個行向量，第i個元素是A的第i列的算術(shù)平均值 median(A)：返回一個行向量，第i個元素是A的第i列的中值,4. 累加和與累乘積,求解元素的累加和與累乘積的函數(shù)為cumsum和cumprod函數(shù)： cumsum(A)：返回一個矩陣，第i列是A的第i列的累加和向量 cumprod(A)：返回一個矩陣，第i列是A的第i列的累乘積向量,見P58 例3-19,5. 排序與標(biāo)準(zhǔn)方差,(1)求標(biāo)準(zhǔn)方差

11、 計算數(shù)據(jù)序列的標(biāo)準(zhǔn)方差的函數(shù)std，std(A)返回一個行向量，它的各個元素便是矩陣A各列的標(biāo)準(zhǔn)方差 Y=std(A) 求各列元素的標(biāo)準(zhǔn)方差,(2). 排序 排序函數(shù)是sort，sort(A)函數(shù)對矩陣A的各列進行排序： Y=sort(A) 對A的各列進行排序,見P58 例3-19,DX=diff(X) 計算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,n-1。,DX=diff(X,n) 計算向量X的n階向前差分，如 diff(X,2)=diff(diff(X)。,DX=diff(A,n,dim) 計算矩陣A的n階差分，dim=1時(缺省狀態(tài))，按列計算差分；dim=2，

12、按行計算差分。,6. 差分,eg3-9,四、插值和擬合,插值的定義是對某些集合給定的數(shù)據(jù)點之間函數(shù)的估值方法。 當(dāng)不能很快地求出所需中間點的函數(shù)時，插值是一個非常有價值的工具。 特點： 利用已知點確定未知點 粗糙 精確 集合大的 簡化的 matlab提供了一維、二維、 三次樣條等許多插值選擇,插值函數(shù)interp1，其調(diào)用格式為： Yi=interp1(X,Y,Xi,method) 函數(shù)根據(jù)X,Y的值，計算函數(shù)在Xi處的值。X,Y是兩個等長的已知向量，分別描述采樣點和樣本值，Xi是一個向量或標(biāo)量，描述欲插值的點，Yi是一個與Xi等長的插值結(jié)果。,1、一維插值函數(shù),method是插值方法，可用的

13、方法有： linear 默認(rèn)方法，線性插值 nearest 鄰近點插值 spline 三次樣條插值 cubic 三次插值，要求x的值等距離, MATLAB中有一個專門的3次樣條插值函數(shù)Yi=spline(X,Y,Xi)，其功能及使用方法與函數(shù)Yi=interp1(X,Y,Xi,spline)基本相同。,注意：1.所有插值方法均要求x是單調(diào)的 2.Xi的取值范圍不能超出X的范圍，否則會給出“NaN”錯誤,eg3-10,優(yōu)缺點 見P63,例3-11 某觀測站測得某日6:00時至18:00時之間每隔2小時的室內(nèi)外溫度()，用3次樣條插值分別求得該日室內(nèi)外6:30至17:30時之間每隔2小時各點的近似

14、溫度(),設(shè)時間變量h為一行向量，溫度變量t為一個兩列矩陣，其中第一列存放室內(nèi)溫度，第二列儲存室外溫度 命令如下： h =6:2:18; t=18,20,22,25,30,28,24;15,19,24,28,34,32,30; XI =6.5:2:17.5 YI=interp1(h,t,XI,spline) %用3次樣條插值計算,eg3-11,2、二維數(shù)據(jù)插值,二維插值函數(shù)interp2，其調(diào)用格式為： Zi=interp2(X,Y,Z,Xi,Yi,method),例3-12 某實驗對一根長10米的鋼軌進行熱源的溫度傳播測試，用x表示測量點0:2.5:10(米)，用h表示測量時間0:30:60

15、(秒)，用T表示測試所得各點的溫度() ，用線性插值求出在一分鐘內(nèi)每隔20秒、鋼軌每隔1米處的溫度TI。 命令如下： x=0:2.5:10; h=0:30:60; T=95,14,0,0,0;88,48,32,12,6;67,64,54,48,41; xi=0:10; hi=0:20:60; TI=interp2(x,h,T,xi,hi),eg3-12,另見P64 例3-25,3、擬合,polyfit函數(shù)可用來求得最小二乘擬合多項式的系數(shù)，其調(diào)用格式為： P=polyfit(x,y,n) 函數(shù)根據(jù)采樣點X和采樣點函數(shù)值Y，產(chǎn)生一個n次多項式P,其中x,y是兩個等長的向量，P是一個長度為n+1的

16、多項式系數(shù)向量。 可用polyval函數(shù)按所得多項式計算給出點上的函數(shù)近似值,例3-13 求所給樣本值的3階擬合多項式 x0=0:0.1:1； y0=-.447 1.978 3.11 5.25 5.02 4.66 4.01 4.58 3.45 5.35 9.22; p=polyfit(x0,y0,3) p = 56.6915 -87.1174 40.0070 -0.9043 xx=0:0.01:1;yy=polyval(p,xx); plot(xx,yy,-b,x0,y0,or),eg3-13,小 結(jié) 本節(jié)介紹了matlab語言的數(shù)值運算 功能，通過學(xué)習(xí)應(yīng)該掌握： 多項式運算 線性方程組 數(shù)據(jù)

17、分析 擬合和插值,作 業(yè) P70 3.1 3.2 3.4 3.5 3.73.12 3.15,五、快速傅里葉變換,一維離散傅立葉變換函數(shù)，其調(diào)用格式與功能為： (1) fft(X)：返回向量X的離散傅立葉變換。設(shè)X的長度(即元素個數(shù))為N，若N為2的冪次，則為以2為基數(shù)的快速傅立葉變換，否則為運算速度很慢的非2冪次的算法。對于矩陣X，fft(X)應(yīng)用于矩陣的每一列。 (2) fft(X,N)：計算N點離散傅立葉變換。它限定向量的長度為N，若X的長度小于N，則不足部分補上零；若大于N，則刪去超出N的那些元素。對于矩陣X，它同樣應(yīng)用于矩陣的每一列，只是限定了向量的長度為N。 (3) fft(X,di

18、m)或fft(X,N,dim)：這是對于矩陣而言的函數(shù)調(diào)用格式，前者的功能與fft(X)基本相同，而后者則與fft(X,N)基本相同。只是當(dāng)參數(shù)dim=1時，該函數(shù)作用于X的每一列；當(dāng)dim=2時，則作用于X的每一行。,五、快速傅里葉變換,值得一提的是，當(dāng)已知給出的樣本數(shù)N0不是2的冪次時，可以取一個N使它大于N0且是2的冪次，然后利用函數(shù)格式fft(X,N)或fft(X,N,dim)便可進行快速傅立葉變換。這樣，計算速度將大大加快。 相應(yīng)地，一維離散傅立葉逆變換函數(shù)是ifft。ifft(F)返回F的一維離散傅立葉逆變換；ifft(F,N)為N點逆變換；ifft(F,dim)或ifft(F,N,dim)則由N或dim確定逆變換的點數(shù)或操作方向。,例3-14 給定數(shù)學(xué)函數(shù) x(t)=12sin(210t+/4)+5cos(240t) 取N=128，試對t從01秒采樣，用fft作快速傅立葉變換，繪制相應(yīng)的振幅-頻率圖 在01秒時間范圍內(nèi)采樣128點，從而可以確定采樣周期和采樣頻率。由于離散