1、第六章 MATLAB的數(shù)值計(jì)算, matlab 具有出色的數(shù)值計(jì)算能力，占據(jù)世界上數(shù)值計(jì)算軟件的主導(dǎo)地位,一、代數(shù)方程組求解,matlab中有兩種除運(yùn)算左除和右除。 對(duì)于方程ax=b，a 為amn矩陣，有三種情 況： 當(dāng)n=m時(shí)，此方程成為“恰定”方程 當(dāng)nm時(shí)，此方程成為“欠定”方程 matlab定義的除運(yùn)算可以很方便地解上 述三種方程,1.恰定方程組的解,方程ax=b(a為非奇異) x=a-1 b 矩陣逆 兩種解: x=inv(a)b 采用求逆運(yùn)算解方程 x=ab 采用左除運(yùn)算解方程,方程ax=b a=1 2;2 3;b=8;13; x=inv(a)*b x=ab x = x = 2.00

2、 2.00 3.00 3.00,=,a x = b,例: x1+2x2=8 2x1+3x2=13,2.超定方程組的解,方程 ax=b ,mn時(shí)此時(shí)不存在唯一解。 方程解 (a a)x=a b x=(a a)-1 a b 求逆法 x=ab matlab用最小二乘法找一 個(gè)準(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,=,a x = b,3.欠定方程組的解,當(dāng)方程數(shù)少于未知量個(gè)數(shù)時(shí),即不定 情況,有無(wú)窮多個(gè)解存在。 ma

3、tlab可求出兩個(gè)解： 一個(gè)是用除法求的解，是具有最多零元素的解 一個(gè)是具有最小長(zhǎng)度或范數(shù)的解，這個(gè)解是基于偽逆pinv求得的。,x1+2x2+3x3=1 2x1+3x2+4x3=2 a=1 2 3;2 3 4;b=1;2; x=ab x=pinv(a)b x = x = 1.00 0.83 0 0.33 0 -0.17,a x = b,借助matlab提供的函數(shù)，處理多項(xiàng)式是一件非常簡(jiǎn)單的事情，很容易對(duì)多項(xiàng)式進(jìn)行積分、微分以及求根的操作。,一元多項(xiàng)式在代數(shù)中占有非常重要的地位。在實(shí)際應(yīng)用中如對(duì)實(shí)驗(yàn)數(shù)據(jù)的插值、微商和曲線擬合等，都要大量用到多項(xiàng)式；在矩陣分析時(shí)，也要用到一元多項(xiàng)式的概念。 多項(xiàng)

4、式函數(shù)是形式最簡(jiǎn)單的函數(shù)，也是最容易計(jì)算的函數(shù)，從理論上講，它可以表示絕大多數(shù)復(fù)雜函數(shù)。 在許多計(jì)算機(jī)的計(jì)算和編程中，很多函數(shù)值如sin(x),cos(x)等的計(jì)算都是先將函數(shù)進(jìn)行Tailor展開為多項(xiàng)式進(jìn)行逼近計(jì)算的，并且都能達(dá)到很高的精度。,二、 多項(xiàng)式運(yùn)算及其求根,鑒于MATLAB無(wú)零下標(biāo)，故把多項(xiàng)式的一般形式表達(dá)為：,在MATLAB中可以用長(zhǎng)度n+1的行向量表示為：,1. 多項(xiàng)式求根,命令格式：x=roots(A)。這里A為多項(xiàng)式的系數(shù)A(1),A(2),A(N),A(N+1)；解得的根賦值給數(shù)組X，即X(1),X(2), ,X(N)。 【例6】試用ROOTS函數(shù)求多項(xiàng)式x4+8x3-

5、10的根 這是一個(gè)4次多項(xiàng)式，它的五個(gè)系數(shù)依次為：1,8,0,0,-10。下面先產(chǎn)生多項(xiàng)式系數(shù)的向量A，然后求根： A=1 8 0 0 -10 A = 1 8 0 0 -10 x=roots(A),x = -8.0194 -0.5075 + 0.9736i -0.5075 - 0.9736i 1.0344,2. 多項(xiàng)式的建立,若已知多項(xiàng)式的全部根，則可以用POLY函數(shù)建立起該多項(xiàng)式；也可以用POLY函數(shù)求矩陣的特征多項(xiàng)式。POLY函數(shù)是一個(gè)MATLAB程序，調(diào)用它的命令格式是： A=poly(x) 若x為具有N個(gè)元素的向量，則poly(x)建立以x為其根的多項(xiàng)式，且將該多項(xiàng)式的系數(shù)賦值給向量A

6、。在此種情況下，POLY與ROOTS互為逆函數(shù)；若x為NN的矩陣x，則poly(x)返回一個(gè)向量賦值給A，該向量的元素為矩陣x的特征多項(xiàng)式之系數(shù)：A(1),A(2),A(N),A(N+1)。,【例7】試用POLY函數(shù)對(duì)例7.8所求得的根，建立相應(yīng)的多項(xiàng)式。 x=-8.0194 -0.5075 + 0.9736i -0.5075 - 0.9736i 1.0344; z=poly(x) z = 1.0000 8.0000 0.0000 0.0000 -9.9996,3. 求多項(xiàng)式的值,POLYVAL函數(shù)用來(lái)求代數(shù)多項(xiàng)式的值，調(diào)用的命令格式為： Y=polyval(A,x) 本命令將POLYVAL函

7、數(shù)返回的多項(xiàng)式的值賦值給Y。若x為一數(shù)值，則Y也為一數(shù)值；若x為向量或矩陣，則對(duì)向量或矩陣中的每個(gè)元素求其多項(xiàng)式的值。,【例8】以4次多項(xiàng)式、分別取x=1.2和下面的矩陣的23個(gè)元素為自變量計(jì)算該多項(xiàng)式的值。 A=1 8 0 0 -10; % 例7.8的4次多項(xiàng)式系數(shù) x=1.2; % 取自變量為一數(shù)值 y1=polyval(A,x) y1 = -5.8976 x=-1 1.2 -1.4;2 -1.8 1.6 % 給出一個(gè)矩陣x y1=polyval(A,x) y1 = -17.0000 5.8976 -28.1104 70.0000 -46.1584 29.3216,4. 多項(xiàng)式的四則運(yùn)算,

8、(1)多項(xiàng)式加、減 對(duì)于次數(shù)相同的若干個(gè)多項(xiàng)式，可直接對(duì)多項(xiàng)式系數(shù)向量進(jìn)行加、減的運(yùn)算。如果多項(xiàng)式的次數(shù)不同，則應(yīng)該把低次的多項(xiàng)式系數(shù)不足的高次項(xiàng)用零補(bǔ)足，使同式中的各多項(xiàng)式具有相同的次數(shù)。,(2)多項(xiàng)式乘法 若A、B是由多項(xiàng)式系數(shù)組成的向量，則CONV函數(shù)將返回這兩個(gè)多項(xiàng)式的乘積。調(diào)用它的命令格式為： C=conv(A,B) 命令的結(jié)果C為一個(gè)向量，由它構(gòu)成一個(gè)多項(xiàng)式。,求4次多項(xiàng)式與多項(xiàng)式2x2-x+3的乘積。 A=1 8 0 0 -10; B=2 -1 3 B = 2 -1 3 C=conv(A,B) C = 2 15 -5 24 -20 10 -30 本例的運(yùn)行結(jié)果是求得一個(gè)6次多項(xiàng)式

9、 2x6+15x5-5x4+24x3-20 x2+10 x-30,(3)多項(xiàng)式除法 當(dāng)A、B是由多項(xiàng)式系數(shù)組成的向量時(shí)，DECONV函數(shù)用來(lái)對(duì)兩個(gè)多項(xiàng)式作除法運(yùn)算。調(diào)用的命令格式為： Q,r=deconv(A,B) 本命令的結(jié)果：多項(xiàng)式A除以多項(xiàng)式B獲商多項(xiàng)式賦予Q（也為多項(xiàng)式系數(shù)向量）；獲余項(xiàng)多項(xiàng)式賦予r（其系數(shù)向量的長(zhǎng)度與被除多項(xiàng)式相同，通常高次項(xiàng)的系數(shù)為0）。 DECONV是CONV的逆函數(shù)，即有A=conv(B,Q)+r。,【例10】試用4次多項(xiàng)式與多項(xiàng)式2x2-x+3相除。 A=1 8 0 0 -10; B=2 -1 3; P,r=deconv(A,B) P = 0.5000 4.2

10、500 1.3750 r = 0 0 0 -11.3750 -14.1250 商多項(xiàng)式P為0.5x2+4.25x+1.375， 余項(xiàng)多項(xiàng)式r為 -11.375x-14.125。,polyder( )多項(xiàng)式微分,k = polyder(p)求多項(xiàng)式p的微分多項(xiàng)式,k = polyder(a,b)求a x b的微分,q,d = polyder(b,a)：求有理分式b/a的微分 q：分子多項(xiàng)式 d：分母多項(xiàng)式,(4).多項(xiàng)式微分,例：a=1 2 3 4 5; poly2str(a,x) ans = x4 + 2 x3 + 3 x2 + 4 x + 5 b=polyder(a) b = 4 6 6 4

11、 poly2str(b,x) ans =4 x3 + 6 x2 + 6 x + 4,10、多項(xiàng)式的積分,g = 1,6,20,48,69,72,44 h = polyder(g) q = polyint(h,44),q = polyint(p,k)對(duì)多項(xiàng)式p積分，k為積分常數(shù) q = polyint(p)對(duì)多項(xiàng)式p積分，積分常數(shù)為0,也可以用p./(length(p):-1:1),k來(lái)完成積分 其中k為常數(shù)項(xiàng) 例：a=1 2 3 4 5; b=polyder(a) b = 4 6 6 4 c= b./(length(b):-1:1),0 c = 1 2 3 4 0,三、 擬合與插值,1. 多項(xiàng)

12、式擬合：用一個(gè)多項(xiàng)式來(lái)逼近一組給定的數(shù)據(jù)。擬合的準(zhǔn)則是最小二乘法 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),polyfit函數(shù)來(lái)求得最小二乘擬合多項(xiàng)式的系數(shù)， polyval函數(shù)按所得的多項(xiàng)式計(jì)算所給出的點(diǎn)上的函數(shù)近似值。 polyfit函數(shù)的調(diào)用格式為： P,S=polyf

13、it(X,Y,m) 函數(shù)根據(jù)采樣點(diǎn)X和采樣點(diǎn)函數(shù)值Y，產(chǎn)生一個(gè)m次多項(xiàng)式P及其在采樣點(diǎn)的誤差向量S。 其中X,Y是兩個(gè)等長(zhǎng)的向量，P是一個(gè)長(zhǎng)度為m+1的向量，P的元素為多項(xiàng)式系數(shù) polyval(P,x) 按多項(xiàng)式的系數(shù)計(jì)算x點(diǎn)多項(xiàng)式的值.,已知數(shù)據(jù)表t,y， t=0,1.7,2.9,3.5,4.8,6.9,7.3,8.6,9.4,10.2 y=5.6,7.8,14,18,19.6,18.6,17.6,15.4,12.5,14 試求2次擬合多項(xiàng)式p(t)，然后求ti=1,1.5,2,2.5,9.5,10各點(diǎn)的函數(shù)近似值。 解：源程序如下 T=0,1.7,2.9,3.5,4.8,6.9,7.3,

14、8.6,9.4,10.2 Y=5.6,7.8,14,18,19.6,18.6,17.6,15.4,12.5,14 P=polyfit(T,Y,2) ti=1:0.5:10 yi=polyval(P,ti) plot(T,Y,r,ti,yi,b),多項(xiàng)式擬合時(shí)的注意事項(xiàng),在進(jìn)行曲線擬合時(shí)對(duì)多項(xiàng)式階次的選擇是任意的。,雖然高階的多項(xiàng)式可以更準(zhǔn)確地?cái)M合數(shù)據(jù)（僅指該曲線與給定的數(shù)據(jù)之間的均方誤差最?。谶M(jìn)行曲線擬合時(shí)，并不需要采用太高階的多項(xiàng)式，這主要基于以下原因：,1、越是高階的多項(xiàng)式其數(shù)值特性越差，計(jì)算起來(lái)也越耗時(shí)；,2、隨著多項(xiàng)式階次的升高，擬合的曲線變的越來(lái)越不平滑，通常會(huì)出現(xiàn)用戶不愿意看

15、到的局部波形；,3、由于數(shù)據(jù)本身的近似性，因此在進(jìn)行數(shù)據(jù)擬合時(shí)沒有必要僅僅考慮使擬合的曲線無(wú)限接近數(shù)據(jù)點(diǎn)，而要在曲線的階次合均方誤差之間綜合考慮，因?yàn)樵绞歉唠A的多項(xiàng)式在物理實(shí)現(xiàn)時(shí)越困難。,從數(shù)學(xué)原理上，n+1個(gè)數(shù)據(jù)點(diǎn)可以惟一定義一個(gè)n次曲線（或n階多項(xiàng)式）。,2.插值 插值是根據(jù)數(shù)據(jù)點(diǎn)的規(guī)律，找到一個(gè)多項(xiàng)式可以連接兩個(gè)點(diǎn)，插值得出兩個(gè)數(shù)據(jù)點(diǎn)之間的數(shù)值。 當(dāng)不能很快地求出所需中間點(diǎn)的函數(shù)時(shí)，插值是一個(gè)非常有價(jià)值的工具。 Matlab提供了一維、二維、 三次樣條等許多插值選擇,y1=interp1(x,y,x1,方法); 例： x = 0:10; y = sin(x); xx = 0:.25:10

16、; yy = spline(x,y,xx); yy=interp1(x,y,xx,spline); plot(x,y,o,xx,yy),Y1=interp1(X,Y,X1,method) 函數(shù)根據(jù)X,Y的值，計(jì)算函數(shù)在X1處的值。 X,Y是兩個(gè)等長(zhǎng)的已知向量，分別描述采樣點(diǎn)和樣本值， X1是一個(gè)向量或標(biāo)量，描述欲插值的點(diǎn)， Y1是一個(gè)與X1等長(zhǎng)的插值結(jié)果。 method是插值方法，允許的取值有 linear、nearest、cubic、spline。 注意：X1的取值范圍不能超出X的給定范圍，否則，會(huì)給出“NaN”錯(cuò)誤。,用不同的插值方法計(jì)算在/2點(diǎn)的值。 x=0,0.5,0.8,1.5,1.

17、9,2.3 y=4,5,8,9,6,3 y1=spline(x,y,pi/2) y1 = 8.5931e+000 y2=interp1(x,y,pi/2,linear) y2 = 8.4690e+000,某觀測(cè)站測(cè)得某日6:00時(shí)至18:00時(shí)之間每隔2小時(shí)的室內(nèi)外溫度()，用3次樣條插值分別求得該日室內(nèi)外6:30至17:30時(shí)之間每隔2小時(shí)各點(diǎn)的近似溫度()。 解：設(shè)時(shí)間變量h為一行向量，溫度變量t為一個(gè)兩列矩陣，其中第一列存放室內(nèi)溫度，第二列儲(chǔ)存室外溫度。命令如下： h =6:2:18; t=18,20,22,25,30,28,24;15,19,24,28,34,32,30; XI =6.

18、5:2:17.5 YI=interp1(h,t,XI,spline) %用3次樣條插值計(jì)算,二維數(shù)據(jù)插值 在MATLAB中，提供了解決二維插值問題的函數(shù)interp2，其調(diào)用格式為： Z1=interp2(X,Y,Z,X1,Y1,method) X,Y是兩個(gè)向量，分別描述兩個(gè)參數(shù)的采樣點(diǎn)， Z是與參數(shù)采樣點(diǎn)對(duì)應(yīng)的函數(shù)值， X1,Y1是兩個(gè)向量或標(biāo)量，描述欲插值的點(diǎn)。 Z1是根據(jù)相應(yīng)的插值方法得到的插值結(jié)果。 method的取值 linear、nearest、cubic、spline X,Y,Z也可以是矩陣形式。 同樣，X1,Y1的取值范圍不能超出X,Y的給定范圍，否則，會(huì)給出“NaN”錯(cuò)誤。,

19、某實(shí)驗(yàn)對(duì)一根長(zhǎng)10米的鋼軌進(jìn)行熱源的溫度傳播測(cè)試。用x表示測(cè)量點(diǎn)0:2.5:10(米)，用h表示測(cè)量時(shí)間0:30:60(秒)，用T表示測(cè)試所得各點(diǎn)的溫度()。試用線性插值求出在一分鐘內(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),四、泛函命令,函數(shù)句柄 方便地實(shí)現(xiàn)函數(shù)間互相調(diào)用 兼容函數(shù)加載的所有方式 拓寬子函數(shù)包括局部函數(shù)的使用范圍 提高函數(shù)調(diào)用的可靠性 減少程序設(shè)計(jì)中的冗余

20、提高重復(fù)執(zhí)行的效率 數(shù)組、結(jié)構(gòu)數(shù)組、細(xì)胞型數(shù)組結(jié)合定義數(shù)據(jù),函數(shù)句柄的創(chuàng)建和顯示 定義函數(shù)句柄的定義只需在提示符后添加相應(yīng)函數(shù)的函數(shù)名 f_h=plot 函數(shù)句柄的內(nèi)容通過functions顯示 functions(f_h) ans = function: plot type: overloaded file: MATLAB built-in function methods: 1x1 struct,matlab的函數(shù)類型 simple 未加載的內(nèi)部函數(shù)、m文件 overloaded 加載的內(nèi)部函數(shù)、m文件 subfunction 子函數(shù) constructor matlab類的創(chuàng)建函數(shù) pr

21、ivate 局部函數(shù),函數(shù)句柄與函數(shù)名字符串轉(zhuǎn)換 func2str 函數(shù)名字符串轉(zhuǎn)換函數(shù)句柄 funh=str2func(det) funh = det str2func 函數(shù)句柄轉(zhuǎn)換函數(shù)名字符串 func2str(funh) ans =det,(1)內(nèi)聯(lián)函數(shù) inline(string,arg1,arg2) 例：f=inline(sin(x),x);f(1:10) (2)函數(shù)句柄 h=fun 或 h=str2func(fun) 例：f=aaa;f(0:pi/100;2*pi),1、尋優(yōu)函數(shù)求單變量函數(shù)最小值 x,y=fminbnd(h_fun,x1,x2) x,y=fminbnd(fun,x1,x2) 其中：y是計(jì)算出的函數(shù)最小值 x是函數(shù)最小值的自變量值 h_fun是函數(shù)句柄 fun為函數(shù)名字符串 x1、x2是自變量取值范圍,例：x = fminbnd(fun_name,2,5) x = fminbnd(sin(x)+3,2,5) x = fminbnd(x) sin(x)+3,2,5) f=inline(sin(x)+3,x); x = fminbnd(f,2,5),2、尋優(yōu)函數(shù)求多變量函數(shù)無(wú)約束最小值（單純形法） x,y=fminsearch(h_fun,x0) x,y=fminsearch(fun,x0) 其中：y是