數(shù)學(xué)建模之?dāng)?shù)值計(jì)算與數(shù)據(jù)分析

1、 第2章數(shù)值計(jì)算與數(shù)據(jù)分析2.1 基本數(shù)學(xué)函數(shù)函數(shù)sin、sinh功能正弦函數(shù)與雙曲正弦函數(shù)格式Y(jié) = sin(X %計(jì)算參量X(可以是向量、矩陣,元素可以是復(fù)數(shù)中每一個(gè)角度分量的正弦值Y,所有分量的角度單位為弧度。Y = sinh(X %計(jì)算參量X的雙曲正弦值Y注意:sin(pi并不是零,而是與浮點(diǎn)精度有關(guān)的無窮小量eps,因?yàn)閜i僅僅是精確值浮點(diǎn)近似的表示值而已;對于復(fù)數(shù)Z= x+iy,函數(shù)的定義為:sin(x+iy = sin(x*cos(y +i*cos(x*sin(y,2eezsin(iz iz-=,2eezsin(z z-=例2-1x = -pi:0.01:pi; plot(x,s

2、in(xx = -5:0.01:5; plot(x,sinh(x圖形結(jié)果為圖2-1。 圖2-1 正弦函數(shù)與雙曲正弦函數(shù)圖函數(shù)asin、asinh功能反正弦函數(shù)與反雙曲正弦函數(shù)格式Y(jié) = asin(X %返回參量X(可以是向量、矩陣中每一個(gè)元素的反正弦函數(shù)值Y。若X中有的分量處于-1,1之間,則Y = asin(X對應(yīng)的分量處于-/2,/2之間,若X中有分量在區(qū)間-1,1之外,則Y= asin(X對應(yīng)的分量為復(fù)數(shù)。Y = asinh(X %返回參量X中每一個(gè)元素的反雙曲正弦函數(shù)值Y說明 反正弦函數(shù)與反雙曲正弦函數(shù)的定義為:z 1z i ln(i z sin a 2-+-=,z 1z ln(z s

3、inh a 2+=例2-2x = -1:.01:1; plot(x,asin(xx = -5:.01:5; plot(x,asinh(x圖形結(jié)果為圖2-2。 圖2-2 反正弦函數(shù)與反雙曲正弦函數(shù)圖 函數(shù) cos 、cosh功能 余弦函數(shù)與雙曲余弦函數(shù)格式 Y = cos(X %計(jì)算參量X(可以是向量、矩陣,元素可以是復(fù)數(shù)中每一個(gè)角度分量的余弦值Y ,所有角度分量的單位為弧度。我們要指出的是,cos(pi/2并不是精確的零,而是與浮點(diǎn)精度有關(guān)的無窮小量eps ,因?yàn)閜i 僅僅是精確值浮點(diǎn)近似的表示值而已。Y = sinh(X %計(jì)算參量X 的雙曲余弦值Y說明 若X 為復(fù)數(shù)z= x+iy ,則函數(shù)

4、定義為:cos(x+iy = cos(x*cos(y + i*sin(x*sin(y,2e e z cos iz iz -+=,2e e z cosh zz -+= 例2-3x = -pi:0.01:pi; plot(x,cos(xx = -5:0.01:5; plot(x,cosh(x圖形結(jié)果為圖2-3。 圖2-3 余弦函數(shù)與雙曲余弦函數(shù)圖 函數(shù) acos 、acosh功能 反余弦函數(shù)與反雙曲余弦函數(shù) 格式 Y = acos(X %返回參量X(可以是向量、矩陣中每一個(gè)元素的反余弦函數(shù)值Y 。若X 中有的分量處于-1,1之間,則Y = acos(X對應(yīng)的分量處于0,之間,若X 中有分量在區(qū)間-

5、1,1之外,則Y = acos(X對應(yīng)的分量為復(fù)數(shù)。Y = asinh(X %返回參量X 中每一個(gè)元素的反雙曲余弦函數(shù)Y說明 反余弦函數(shù)與反雙曲余弦函數(shù)定義為:z 1i z i ln(i z co s a 2-+-=,1z z ln(z cosh a 2-+=例2-4x = -1:.01:1; plot(x,acos(xx = -5:.01:5; plot(x,acosh(x圖形結(jié)果為圖2-4。 圖2-4 反余弦函數(shù)與反雙曲余弦函數(shù)圖 函數(shù) tan 、tanh功能 正切函數(shù)與雙曲正切函數(shù)格式 Y = tan(X %計(jì)算參量X(可以是向量、矩陣,元素可以是復(fù)數(shù)中每一個(gè)角度分量的正切值Y ,所有角

6、度分量的單位為弧度。我們要指出的是,tan(pi/2并不是精確的零,而是與浮點(diǎn)精度有關(guān)的無窮小量eps ,因?yàn)閜i 僅僅是精確值浮點(diǎn)近似的表示值而已。Y = tanh(X %返回參量X 中每一個(gè)元素的雙曲正切函數(shù)值Y例2-5x = (-pi/2+0.01:0.01:(pi/2-0.01; % 稍微縮小定義域plot(x,tan(xx = -5:0.01:5; plot(x,tanh(x圖形結(jié)果為圖2-5。 圖2-5 正切函數(shù)與雙曲正切函數(shù)圖 函數(shù) atan 、atanh功能 反正切函數(shù)與反雙曲正切函數(shù)格式 Y = atan(X %返回參量X(可以是向量、矩陣中每一個(gè)元素的反正切函數(shù)值Y 。若X

7、 中有的分量為實(shí)數(shù),則Y = atan(X對應(yīng)的分量處于-/2,/2之間。Y = atanh(X %返回參量X 中每一個(gè)元素的反雙曲正切函數(shù)值Y 。說明 反正切函數(shù)與反雙曲正切函數(shù)定義為:z i z i ln 2i z tan a -+=,z1z 1ln 21z tanh a -+= 例2-6x = -20:0.01:20; plot(x,atan(xx = -0.99:0.01:0.99; plot(x,atanh(x 圖形結(jié)果為圖2-6。 圖2-6 反正切函數(shù)與反雙曲正切函數(shù)圖 函數(shù) cot 、coth功能 余切函數(shù)與雙曲余切函數(shù)格式 Y = cot(X %計(jì)算參量X(可以是向量、矩陣,元

8、素可以是復(fù)數(shù)中每一個(gè)角度分量的余切值Y ,所有角度分量的單位為弧度。Y = coth(X %返回參量X 中每一個(gè)元素的雙曲余切函數(shù)值Y例2-7x1 = -pi+0.01:0.01:-0.01; % 去掉奇點(diǎn)x = 0x2 = 0.01:0.01:pi-0.01; % 做法同上plot(x1,cot(x1,x2,cot(x2plot(x1,coth(x1,x2,coth(x2圖形結(jié)果為圖2-7。 圖2-7 余切函數(shù)與雙曲余切函數(shù)圖 函數(shù)acot、acoth功能反余切函數(shù)與反雙曲余切函數(shù)格式Y(jié) = acot(X %返回參量X(可以是向量、矩陣中每一個(gè)元素的反余切函數(shù)Y Y = acoth(X %返

9、回參量X中每一個(gè)元素的反雙曲余切函數(shù)值Y例2-8x1 = -2*pi:pi/30:-0.1; x2 = 0.1:pi/30:2*pi; % 去掉奇異點(diǎn)x = 0plot(x1,acot(x1,x2,acot(x2x1 = -30:0.1:-1.1; x2 = 1.1:0.1:30;plot(x1,acoth(x1,x2,acoth(x2圖形結(jié)果為圖2-8。 圖2-8 反余切函數(shù)與反雙曲余切函數(shù)圖函數(shù)sec、sech功能正割函數(shù)與雙曲正割函數(shù)格式Y(jié) = sec(X %計(jì)算參量X(可以是向量、矩陣,元素可以是復(fù)數(shù)中每一個(gè)角度分量的正割函數(shù)值Y,所有角度分量的單位為弧度。我們要指出的是,sec(pi

10、/2并不是無窮大,而是與浮點(diǎn)精度有關(guān)的無窮小量eps的倒數(shù),因?yàn)閜i僅僅是精確值浮點(diǎn)近似的表示值而已。Y = sech(X %返回參量X中每一個(gè)元素的雙曲正割函數(shù)值Y例2-9x1 = -pi/2+0.01:0.01:pi/2-0.01; % 去掉奇異點(diǎn)x = pi/2x2 = pi/2+0.01:0.01:(3*pi/2-0.01;plot(x1,sec(x1,x2,sec(x2x = -2*pi:0.01:2*pi;plot(x,sech(x圖形結(jié)果為圖2-9。 圖2-9 正割函數(shù)與雙曲正割函數(shù)圖 函數(shù)asec、asech功能反正割函數(shù)與反雙曲正割函數(shù)格式Y(jié) = asec(X %返回參量X(

11、可以是向量、矩陣中每一個(gè)元素的反正割函數(shù)值Y Y = asech(X %返回參量X中每一個(gè)元素的反雙曲正割函數(shù)值Y 例2-10x1 = -5:0.01:-1; x2 = 1:0.01:5;plot(x1,asec(x1,x2,asec(x2x = 0.01:0.001:1; plot(x,asech(x圖形結(jié)果為圖2-10。 圖2-10 反正割函數(shù)與反雙曲正割函數(shù)圖函數(shù)csc、csch功能余割函數(shù)與雙曲余割函數(shù)格式Y(jié) = csc(X %計(jì)算參量X(可以是向量、矩陣,元素可以是復(fù)數(shù)中每一個(gè)角度分量的余割函數(shù)值Y,所有角度分量的單位為弧度。Y = csch(X %返回參量X中每一個(gè)元素的雙曲余割函

12、數(shù)值Y例2-11x1 = -pi+0.01:0.01:-0.01; x2 = 0.01:0.01:pi-0.01; % 去掉奇異點(diǎn)x=0plot(x1,csc(x1,x2,csc(x2plot(x1,csch(x1,x2,csch(x2圖形結(jié)果為圖2-11。 圖2-11 余割函數(shù)與雙曲余割函數(shù)圖函數(shù)acsc、acsch功能反余割函數(shù)與反雙曲余割函數(shù)。 格式 Y = asec(X %返回參量X(可以是向量、矩陣中每一個(gè)元素的反余割函數(shù)值YY = asech(X %返回參量X 中每一個(gè)元素的反雙曲余割函數(shù)值Y 例2-12x1 = -10:0.01:-1.01; x2 = 1.01:0.01:10;

13、 % 去掉奇異點(diǎn)x = 1 plot(x1,acsc(x1,x2,acsc(x2x1 = -20:0.01:-1; x2 = 1:0.01:20; plot(x1,acsch(x1,x2,acsch(x2圖形結(jié)果為圖2-12。 圖2-12 反余割函數(shù)與反雙曲余割函數(shù)圖函數(shù) atan2功能 四象限的反正切函數(shù)格式 P = atan2(Y ,X %返回一與參量X 和Y 同型的、與X 和Y 元素的實(shí)數(shù)部分對應(yīng)的、元素對元素的四象限的反正切函數(shù)陣列P ,其中X 和Y 的虛數(shù)部分將忽略。陣列P 中的元素分布在閉區(qū)間-pi,pi上。特定的象限將取決于sign(Y與sign(X。例2-13z=1+2i; r

14、 = abs(z; theta = atan2(imag(z,real(z z = r *exp(i *theta feather(z;hold on t=0:0.1:2*pi;x=1+sqrt(5*cos(t; y=sqrt(5*sin(t; plot(x,y;axis equal; hold off計(jì)算結(jié)果為:theta =1.1071 z =1.0000 +2.0000i圖形結(jié)果為圖2-13。 圖2-13 四象限的反正切函數(shù)圖函數(shù)fix功能朝零方向取整格式 B = fix(A %對A的每一個(gè)元素朝零的方向取整數(shù)部分,返回與A同維的數(shù)組。對于復(fù)數(shù)參量A,則返回一復(fù)數(shù),其分量的實(shí)數(shù)與虛數(shù)部分

15、分別取原復(fù)數(shù)的、朝零方向的整數(shù)部分。例2-14>>A = -1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i;>>B = fix(A計(jì)算結(jié)果為:B =Columns 1 through 4-1.0000 0 3.0000 5.0000Columns 5 through 67.0000 2.0000 + 3.0000i函數(shù)roud功能朝最近的方向取整。格式Y(jié) = round(X %對X的每一個(gè)元素朝最近的方向取整數(shù)部分,返回與X同維的數(shù)組。對于復(fù)數(shù)參量X,則返回一復(fù)數(shù),其分量的實(shí)數(shù)與虛數(shù)部分分別取原復(fù)數(shù)的、朝最近方向的整數(shù)部分。例2-15&

16、gt;>A = -1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i;>>Y = round(A計(jì)算結(jié)果為:Y =Columns 1 through 4-2.0000 0 3.0000 6.0000Columns 5 through 67.0000 2.0000 + 4.0000i函數(shù)floor功能朝負(fù)無窮大方向取整格式 B = floor(A %對A的每一個(gè)元素朝負(fù)無窮大的方向取整數(shù)部分,返回與A同維的數(shù)組。對于復(fù)數(shù)參量A,則返回一復(fù)數(shù),其分量的實(shí)數(shù)與虛 數(shù)部分分別取原復(fù)數(shù)的、朝負(fù)無窮大方向的整數(shù)部分。例2-16>>A = -1.9

17、, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i;>>F = floor(A計(jì)算結(jié)果為:F =Columns 1 through 4-2.0000 -1.0000 3.0000 5.0000Columns 5 through 67.0000 2.0000 + 3.0000i函數(shù)rem功能求作除法后的剩余數(shù)格式R = rem(X,Y %返回結(jié)果X - fix(X./Y.*Y,其中X、Y應(yīng)為正數(shù)。若X、Y為浮點(diǎn)數(shù),由于計(jì)算機(jī)對浮點(diǎn)數(shù)的表示的不精確性,則結(jié)果將可能是不可意料的。fix(X./Y為商數(shù)X./Y朝零方向取的整數(shù)部分。若X與Y為同符號(hào)的,則rem(X,

18、Y返回的結(jié)果與mod(X,Y相同,不然,若X為正數(shù),則rem(-X,Y = mod(-X,Y - Y。該命令返回的結(jié)果在區(qū)間0,sign(X*abs(Y,若Y中有零分量,則相應(yīng)地返回NaN。例2-17>>X = 12 23 34 45;>>Y = 3 7 2 6;>>R = rem(X,Y計(jì)算結(jié)果為:R =0 2 0 3函數(shù)ceil功能朝正無窮大方向取整格式 B = floor(A % 對A的每一個(gè)元素朝正無窮大的方向取整數(shù)部分,返回與A同維的數(shù)組。對于復(fù)數(shù)參量A,則返回一復(fù)數(shù),其分量的實(shí)數(shù)與虛數(shù)部分分別取原復(fù)數(shù)的、朝正無窮大方向的整數(shù)部分。例2-18>

19、;>A = -1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i;>>B = ceil(A計(jì)算結(jié)果為:B =Columns 1 through 4-1.0000 0 4.0000 6.0000Columns 5 through 67.0000 3.0000 + 4.0000i函數(shù)exp功能以e為底數(shù)的指數(shù)函數(shù)格式Y(jié) = exp(X %對參量X的每一分量,求以e為底數(shù)的指數(shù)函數(shù)Y。X中的分量可以為復(fù)數(shù)。對于復(fù)數(shù)分量如,z = x +i*y,則相應(yīng)地計(jì)算:ez =ex*(cos(y + i*sin(y。 例2-19>>A = -1.9,

20、-0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i;>>Y = exp(A計(jì)算結(jié)果為:Y =1.0e+003 *Columns 1 through 40.0001 0.0008 0.0231 0.2704Columns 5 through 61.0966 -0.0099 - 0.0049i函數(shù)expm功能求矩陣的以e為底數(shù)的指數(shù)函數(shù)格式Y(jié) = expm(X %計(jì)算以e為底數(shù)、x的每一個(gè)元素為指數(shù)的指數(shù)函數(shù)值。若矩陣x有小于等于零的特征值,則返回復(fù)數(shù)的結(jié)果。說明該函數(shù)為一內(nèi)建函數(shù),它有三種計(jì)算算法:(1使用文件expm1.m中的用比例法與二次冪算法得到的Pad近似

21、值;(2使用Taylor級(jí)數(shù)近似展開式計(jì)算,這種計(jì)算在文件expm2.m中。但這種一般計(jì)算方法是不可取的,通常計(jì)算是緩慢且不精確的;(3在文件expm3.m中,先是將矩陣對角線化,再把函數(shù)計(jì)算出相應(yīng)的的特征向量,最后轉(zhuǎn)換過來。但當(dāng)輸入的矩陣沒有與矩陣階數(shù)相同的特征向量個(gè)數(shù)時(shí),就會(huì)出現(xiàn)錯(cuò)誤。例2-20>>A=hilb(4;>>Y = expm(A計(jì)算結(jié)果為:Y =3.2506 1.2068 0.8355 0.64171.2068 1.7403 0.5417 0.42880.8355 0.5417 1.4100 0.33180.6417 0.4288 0.3318 1.27

22、29函數(shù)log功能自然對數(shù),即以e為底數(shù)的對數(shù)。格式Y(jié) = log(X %對參量X中的每一個(gè)元素計(jì)算自然對數(shù)。其中X中的元素可以是復(fù)數(shù)與負(fù)數(shù),但由此可能得到意想不到的結(jié)果。若z = x + i*y,則log對復(fù)數(shù)的計(jì)算如下:log (z = log (abs (z + i*atan2(y,x 例2-21 下面的語句可以得到無理數(shù)的近似值:>>Pi = abs(log(-1計(jì)算結(jié)果為:Pi =3.1416函數(shù)log10功能常用對數(shù),即以10為底數(shù)的對數(shù)。格式Y(jié) = log10(X %計(jì)算X中的每一個(gè)元素的常用對數(shù),若X中出現(xiàn)復(fù)數(shù),則可能得到意想不到的結(jié)果。例2-22 >>

23、L1 = log10(realmax % 由此可得特殊變量realmax的近似值>>L2 = log10(eps % 由此可得特殊變量eps的近似值>>M = magic(4;>>L3 = log10(M計(jì)算結(jié)果為:L1 =308.2547L2 =-15.6536L3 =1.2041 0.3010 0.4771 1.11390.6990 1.0414 1.0000 0.90310.9542 0.8451 0.7782 1.07920.6021 1.1461 1.1761 0函數(shù)sort功能把輸入?yún)⒘恐械脑匕磸男〉酱蟮姆较蛑匦屡帕懈袷?B = sort(A

24、%沿著輸入?yún)⒘緼的不同維的方向、從小到大重新排列A中的元素。A可以是字符串的、實(shí)數(shù)的、復(fù)數(shù)的單元數(shù)組。對于A中完全相同的元素,則按它們在A中的先后位置排列在一塊;若A為復(fù)數(shù)的,則按元素幅值的從小到大排列,若有幅值相同的復(fù)數(shù)元素,則再按它們在區(qū)間-,的幅角從小到大排列;若A中有元素為NaN,則將它們排到最后。若A為向量,則返回從小到大的向量,若A為二維矩陣,則按列的方向進(jìn)行排列;若A為多維數(shù)組,sort(A把沿著第一非單元集的元素象向量一樣進(jìn)行處理。B = sort(A,dim %沿著矩陣A(向量的、矩陣的或多維的中指定維數(shù)dim方向重新排列A中的元素。B,INDEX = sort(A, %輸出

25、參量B的結(jié)果如同上面的情形,輸出INDEX是一等于size(A的數(shù)組,它的每一列是與A中列向量的元素相對應(yīng)的置換向量。若A中有重復(fù)出現(xiàn)的相同的值,則返回保存原來相對位置的索引。例2-23>>A = -1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i;>>B1,INDEX = sort(A>>M = magic(4;>>B2 = sort(M計(jì)算結(jié)果為:B1 =Columns 1 through 4-0.2000 -1.9000 3.1416 2.4000 + 3.6000iColumns 5 through 65.

26、6000 7.0000INDEX =2 13 64 5B2 =4 2 3 15 76 89 11 10 12 16 14 15 13函數(shù)abs功能數(shù)值的絕對值與復(fù)數(shù)的幅值格式Y(jié) = abs(X %返回參量X的每一個(gè)分量的絕對值;若X為復(fù)數(shù)的,則返回每一分量的幅值:abs(X = sqrt(real(X.2+imag(X.2。例2-24>>A = -1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i;>>Y = abs(A計(jì)算結(jié)果為:Y =1.9000 0.2000 3.1416 5.6000 7.0000 4.3267函數(shù)conj功能復(fù)數(shù)的

27、共軛值格式ZC = conj(Z %返回參量Z的每一個(gè)分量的共軛復(fù)數(shù):conj(Z = real(Z - i*imag(Z函數(shù)imag功能復(fù)數(shù)的虛數(shù)部分格式Y(jié) = imag(Z %返回輸入?yún)⒘縕的每一個(gè)分量的虛數(shù)部分。例2-25>>imag(2+3i計(jì)算結(jié)果為:ans =3函數(shù)real功能復(fù)數(shù)的實(shí)數(shù)部分。格式Y(jié) = real(Z %返回輸入?yún)⒘縕的每一個(gè)分量的實(shí)數(shù)部分。例2-26>>real(2+3i計(jì)算結(jié)果為:ans =2函數(shù)angle功能復(fù)數(shù)的相角格式P = angle(Z %返回輸入?yún)⒘縕的每一復(fù)數(shù)元素的、單位為弧度的相角,其值在區(qū)間-,上。說明angle(z =

28、 imag (log(z = atan2 (imag(z,real(z例2-27>>Z =1-i, 2+i, 3-i, 4+i;>>1+2i,2-2i,3+2i,4-2i;>>1-3i,2+3i,3-3i,4+3i;>>1+4i,2-4i,3+4i,4-4i;>>P = angle(Z計(jì)算結(jié)果為:P = -0.7854 0.4636 -0.3218 0.24501.1071 -0.7854 0.5880 -0.4636-1.2490 0.9828 -0.7854 0.64351.3258 -1.1071 0.9273 -0.7854函

29、數(shù)complex功能用實(shí)數(shù)與虛數(shù)部分創(chuàng)建復(fù)數(shù)格式 c = complex(a,b %用兩個(gè)實(shí)數(shù)a,b創(chuàng)建復(fù)數(shù)c=a+bi。輸出參量c與a、b同型(同為向量、矩陣、或多維陣列。該命令比下列形式的復(fù)數(shù)輸入更有用:a + i*b 或a + j*b因?yàn)閕和j可能被用做其他的變量(不等于sqrt(-1,或者a和b不是雙精度的。c = complex(a %輸入?yún)⒘縜作為輸出復(fù)數(shù)c的實(shí)部,其虛部為0:c = a+0*i。例2-28>>a = uint8(1;2;3;4;>>b = uint8(4;3;2;1;>>c = complex(a,b計(jì)算結(jié)果為:c =1.000

30、0 + 4.0000i2.0000 +3.0000i3.0000 + 2.0000i4.0000 + 1.0000i函數(shù)mod功能模數(shù)(帶符號(hào)的除法余數(shù)用法M = mod(X,Y %輸入?yún)⒘縓、Y應(yīng)為整數(shù),此時(shí)返回余數(shù)X -Y.*floor(X./Y,若Y0,或者是X。若運(yùn)算數(shù)x與y有相同的符號(hào),則mod(X,Y等于rem(X,Y。總之,對于整數(shù)x,y,有:mod(-x,y =rem(-x,y+y。若輸入為實(shí)數(shù)或復(fù)數(shù),由于浮點(diǎn)數(shù)在計(jì)算機(jī)上的不精確表示,該操作將導(dǎo)致不可預(yù)測的結(jié)果。例2-29>>M1 = mod(13,5>>M2 = mod(1:5,3>>M3

31、 = mod(magic(3,3計(jì)算結(jié)果為:M1 =3M2 =1 2 0 1 2M3 =2 1 00 2 11 0 2函數(shù)nchoosek功能二項(xiàng)式系數(shù)或所有的組合數(shù)。該命令只有對n<15時(shí)有用。函數(shù) C = nchoosek(n,k %參量n,k為非負(fù)整數(shù),返回n! / ( (n-k! k!,即一次從n個(gè)物體中取出k個(gè)的組合數(shù)。C = nchoosek(v,k %參量v為n維向量,返回一矩陣,其行向量的分量為一次C=n! /性從v個(gè)物體中取k個(gè)物體的組合數(shù)。矩陣 C包含nk( (n-k! k!行與k列。例2-30>>C = nchoosek(2:2:10,4計(jì)算結(jié)果為:C

32、=2 4 6 82 4 6 102 4 8 102 6 8 104 6 8 10函數(shù)rand功能生成元素均勻分布于(0,1上的數(shù)值與陣列用法Y = rand(n %返回n*n階的方陣Y,其元素均勻分布于區(qū)間(0,1。若n不是一標(biāo)量,在顯示一出錯(cuò)信息。Y = rand(m,n、Y = rand(m n %返回階數(shù)為m*n的,元素均勻分布于區(qū)間(0,1上矩陣Y。Y = rand(m,n,p,、Y = rand(m n p %生成階數(shù)m*n*p*的,元素服從均勻分布的多維隨機(jī)陣列Y。Y = rand(size(A %生成一與陣列A同型的隨機(jī)均勻陣列Yrand %該命令在每次單獨(dú)使用時(shí),都返回一隨機(jī)數(shù)

33、(服從均勻分布。s = rand('state' %返回一有35元素的列向量s,其中包含均勻分布生成器的當(dāng)前狀態(tài)。該改變生成器的當(dāng)前的狀態(tài),見表2-1。表2-1 >>R1 = rand(4,5>>a = 10; b = 50;>>R2 = a + (b-a * rand(5 % 生成元素均勻分布于(10,50上的矩陣計(jì)算結(jié)果可能為:R1 =0.6655 0.0563 0.2656 0.5371 0.67970.3278 0.4402 0.9293 0.5457 0.61290.6325 0.4412 0.9343 0.9394 0.39400

34、.5395 0.6501 0.5648 0.7084 0.2206R2 =33.6835 19.8216 36.9436 49.6289 46.467918.5164 34.2597 15.3663 31.0549 49.037719.0026 37.1006 33.6046 39.5361 13.933612.4641 12.9804 35.5420 23.2916 46.830428.5238 48.7418 49.0843 13.0512 10.9265函數(shù)randn 功能生成元素服從正態(tài)分布(N(0,1的數(shù)值與陣列格式Y(jié) = randn(n %返回n*n階的方陣Y,其元素服從正態(tài)分布N

35、(0,1。若n不是一標(biāo)量,則顯示一出錯(cuò)信息。Y = randn(m,n、Y = randn(m n %返回階數(shù)為m*n的,元素均勻分布于區(qū)間(0,1上矩陣Y。Y = randn(m,n,p,、Y = randn(m n p %生成階數(shù)m*n*p*的,元素服從正態(tài)分布的多維隨機(jī)陣列Y。Y = randn(size(A %生成一與陣列A同型的隨機(jī)正態(tài)陣列Yrandn %該命令在每次單獨(dú)使用時(shí),都返回一隨機(jī)數(shù)(服從正態(tài)分布。s = randn('state' %返回一有2元素的向量s,其中包含正態(tài)分布生成器的當(dāng)前狀態(tài)。該改變生成器的當(dāng)前狀態(tài),見表2-2。表2-2 例:>>

36、R1 = rand(4,5>>R2 = 0.6 + sqrt(0.1 * randn(5計(jì)算結(jié)果可能為:R1 =0.2778 0.2681 0.5552 0.5167 0.88210.2745 0.3710 0.1916 0.3385 0.58230.9124 0.5129 0.4164 0.2993 0.05500.4125 0.2697 0.1508 0.9370 0.5878R2 =0.4632 0.9766 0.5410 0.6360 0.69310.0733 0.9760 0.8295 0.9373 0.17750.6396 0.5881 0.4140 0.6187 0.

37、82590.6910 0.7035 1.2904 0.5698 1.11340.2375 0.6552 0.5569 0.3368 0.38122.2 插值、擬合與查表插值法是實(shí)用的數(shù)值方法,是函數(shù)逼近的重要方法。在生產(chǎn)和科學(xué)實(shí)驗(yàn)中,自變量x與因變量y的函數(shù)y = f(x的關(guān)系式有時(shí)不能直接寫出表達(dá)式,而只能得到函數(shù)在若干個(gè)點(diǎn)的函數(shù)值或?qū)?shù)值。當(dāng)要求知道觀測點(diǎn)之外的函數(shù)值時(shí),需要估計(jì)函數(shù)值在該點(diǎn)的值。如何根據(jù)觀測點(diǎn)的值,構(gòu)造一個(gè)比較簡單的函數(shù)y=(x,使函數(shù)在觀測點(diǎn)的值等于已知的數(shù)值或?qū)?shù)值。用簡單函數(shù)y=(x在點(diǎn)x處的值來估計(jì)未知函數(shù)y=f(x在x點(diǎn)的值。尋找這樣的函數(shù)(x,辦法是很多的。(

38、x可以是一個(gè)代數(shù)多項(xiàng)式,或是三角多項(xiàng)式,也可以是有理分式;(x可以是任意光滑(任意階導(dǎo)數(shù)連續(xù)的函數(shù)或是分段函數(shù)。函數(shù)類的不同,自然地有不同的逼近效果。在許多應(yīng)用中,通常要用一個(gè)解析函數(shù)(一、二元函數(shù)來描述觀測數(shù)據(jù)。 根據(jù)測量數(shù)據(jù)的類型:1.測量值是準(zhǔn)確的,沒有誤差。2.測量值與真實(shí)值有誤差。這時(shí)對應(yīng)地有兩種處理觀測數(shù)據(jù)方法:1.插值或曲線擬合。2.回歸分析(假定數(shù)據(jù)測量是精確時(shí),一般用插值法,否則用曲線擬合。命令1interp1功能一維數(shù)據(jù)插值(表格查找。該命令對數(shù)據(jù)點(diǎn)之間計(jì)算內(nèi)插值。它找出一元函數(shù)f(x在中間點(diǎn)的數(shù)值。其中函數(shù)f(x由所給數(shù)據(jù)決定。各個(gè)參量之間的關(guān)系示意圖為圖2-14。 圖2

39、-14 數(shù)據(jù)點(diǎn)與插值點(diǎn)關(guān)系示意圖格式y(tǒng)i = interp1(x,Y,xi %返回插值向量yi,每一元素對應(yīng)于參量xi,同時(shí)由向量x與Y的內(nèi)插值決定。參量x指定數(shù)據(jù)Y的點(diǎn)。若Y為一矩陣,則按Y的每列計(jì)算。yi是階數(shù)為length(xi*size(Y,2的輸出矩陣。yi = interp1(Y,xi %假定x=1:N,其中N為向量Y的長度,或者為矩陣Y的行數(shù)。yi = interp1(x,Y,xi,method %用指定的算法計(jì)算插值:nearest:最近鄰點(diǎn)插值,直接完成計(jì)算;linear:線性插值(缺省方式,直接完成計(jì)算;spline:三次樣條函數(shù)插值。對于該方法,命令interp1調(diào)用函數(shù)

40、spline、ppval、mkpp、umkpp。這些命令生成一系列用于分段多項(xiàng)式操作的函數(shù)。命令spline用它們執(zhí)行三次樣條函數(shù)插值;pchip:分段三次Hermite插值。對于該方法,命令interp1調(diào)用函數(shù)pchip,用于對向量x與y執(zhí)行分段三次內(nèi)插值。該方法保留單調(diào)性與數(shù)據(jù)的外形;cubic:與pchip操作相同;v5cubic:在MATLAB 5.0中的三次插值。對于超出x范圍的xi的分量,使用方法nearest、linear、v5cubic的插值算法,相應(yīng)地將返回NaN。對其他的方法,interp1將對超出的分量執(zhí)行外插值算法。 yi = interp1(x,Y,xi,metho

41、d,'extrap' %對于超出x范圍的xi中的分量將執(zhí)行特殊的外插值法extrap。yi = interp1(x,Y,xi,method,extrapval %確定超出x范圍的xi中的分量的外插值extrapval,其值通常取NaN或0。例2-31>>x = 0:10; y = x.*sin(x;>>xx = 0:.25:10; yy = interp1(x,y,xx;>>plot(x,y,'kd',xx,yy插值圖形為圖2-15。例2-32>> year = 1900:10:2010;>> prod

42、uct = 75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 256.344 267.893 ;>>p1995 = interp1(year,product,1995>>x = 1900:1:2010;>>y = interp1(year,product,x,'pchip'>>plot(year,product,'o',x,y插值結(jié)果為:p1995 =252.9885插值圖形為圖2-16。 圖2-15 一元函數(shù)

43、插值圖形圖2-16 離散數(shù)據(jù)的一維插值圖命令2 interp2功能二維數(shù)據(jù)內(nèi)插值(表格查找格式ZI = interp2(X,Y,Z,XI,YI %返回矩陣ZI,其元素包含對應(yīng)于參量XI與YI(可以是向量、或同型矩陣的元素,即Zi(i,jXi(i,j,yi(i,j。用戶可以輸入行向量和列向量Xi與Yi,此時(shí),輸出向量Zi與矩陣meshgrid(xi,yi是同型的。同時(shí)取決于由輸入矩陣X、Y與Z確定的二維函數(shù)Z=f(X,Y。參量X與Y必須是單調(diào)的,且相同的劃分格式,就像由命令meshgrid生成的一樣。若Xi與Yi中有在X與Y范圍之外的點(diǎn),則相應(yīng)地返回nan(Not a Number。ZI = i

44、nterp2(Z,XI,YI %缺省地,X=1:n、Y=1:m,其中m,n=size(Z。再按第一種情形進(jìn)行計(jì)算。ZI = interp2(Z,n %作n 次遞歸計(jì)算,在Z 的每兩個(gè)元素之間插入它們的二維插值,這樣,Z 的階數(shù)將不斷增加。interp2(Z等價(jià)于interp2(z,1。ZI = interp2(X,Y ,Z,XI,YI,method %用指定的算法method 計(jì)算二維插值:linear:雙線性插值算法(缺省算法;nearest :最臨近插值;spline :三次樣條插值;cubic :雙三次插值。例2-33:>>X,Y = meshgrid(-3:.25:3;&g

45、t;>Z = peaks(X,Y; >>XI,YI = meshgrid(-3:.125:3;>>ZZ = interp2(X,Y ,Z,XI,YI;>>surfl(X,Y ,Z;hold on;>>surfl(XI,YI,ZZ+15>>axis(-3 3 -3 3 -5 20;shading flat>>hold off插值圖形為圖2-17。例2-34>>years = 1950:10:1990;>>service = 10:10:30;>>wage = 150.697 199.

46、592 187.625179.323 195.072 250.287203.212 179.092 322.767226.505 153.706 426.730249.633 120.281 598.243;>>w = interp2(service,years,wage,15,1975 插值結(jié)果為:w =190.6288命令3 interp3功能 三維數(shù)據(jù)插值(查表格式 VI = interp3(X,Y ,Z,V ,XI,YI,ZI %找出由參量X,Y,Z 決定的三元函數(shù)V=V(X,Y,Z在點(diǎn)(XI,YI,ZI的值。參量XI,YI,ZI 是同型陣列或向量。若向量參量XI,YI,Z

47、I 是不同長度,不同方向(行或列的向量,這時(shí)輸出參量VI 與Y1,Y2,Y3為同型矩陣。其中Y1,Y2,Y3為用命令meshgrid(XI,YI,ZI生成的同型陣列。若插值點(diǎn)(XI,YI,ZI中有位于點(diǎn)(X,Y ,Z之外的點(diǎn),則相應(yīng)地返回特殊變量值NaN 。VI = interp3(V ,XI,YI,ZI %缺省地,X=1:N ,Y=1:M ,Z=1:P ,其中,M,N,P=size(V,再按上面的情形計(jì)算。VI = interp3(V ,n %作n 次遞歸計(jì)算,在V 的每兩個(gè)元素之間插入它們的三維插值。這樣,V 的階數(shù)將不斷增加。interp3(V等價(jià)于interp3(V,1。圖2-17 二

48、維插值圖 VI = interp3(,method %用指定的算法method作插值計(jì)算:linear:線性插值(缺省算法;cubic:三次插值;spline:三次樣條插值;nearest:最鄰近插值。說明在所有的算法中,都要求X,Y,Z是單調(diào)且有相同的格點(diǎn)形式。當(dāng)X,Y,Z是等距且單調(diào)時(shí),用算法*linear,*cubic,*nearest,可得到快速插值。例2-35>>x,y,z,v = flow(20;>>xx,yy,zz = meshgrid(.1:.25:10, -3:.25:3, -3:.25:3;>>vv = interp3(x,y,z,v,x

49、x,yy,zz;>>slice(xx,yy,zz,vv,6 9.5,1 2,-2 .2; shading interp;colormap cool插值圖形為圖2-18。 圖2-18 三維插值圖命令4 interpft功能用快速Fourier算法作一維插值格式y(tǒng) = interpft(x,n %返回包含周期函數(shù)x在重采樣的n個(gè)等距的點(diǎn)的插值y。若length(x=m,且x有采樣間隔dx,則新的y的采樣間隔dy=dx*m/n。注意的是必須nm。若x為一矩陣,則按x的列進(jìn)行計(jì)算。返回的矩陣y有與x相同的列數(shù),但有n行。y = interpft(x,n,dim %沿著指定的方向dim進(jìn)行計(jì)

50、算命令5 griddata功能數(shù)據(jù)格點(diǎn)格式ZI = griddata(x,y,z,XI,YI %用二元函數(shù)z=f(x,y的曲面擬合有不規(guī)則的數(shù)據(jù)向量x,y,z。griddata將返回曲面z在點(diǎn)(XI,YI處的插值。曲面總是經(jīng)過這些數(shù)據(jù)點(diǎn)(x,y,z的。輸入?yún)⒘?XI,YI通常是規(guī)則的格點(diǎn)(像用命令meshgrid生成的一樣。XI可以是一行向量,這時(shí)XI指定一有常數(shù)列向量的矩陣。類似地,YI可以是一列向量,它指定一有常數(shù)行向量的矩陣。XI,YI,ZI = griddata(x,y,z,xi,yi %返回的矩陣ZI含義同上,同時(shí),返回的矩陣XI,YI是由行向量xi與列向量yi用命令 meshgri

51、d 生成的。 = griddata(,method %用指定的算法method 計(jì)算:linear:基于三角形的線性插值(缺省算法;cubic: 基于三角形的三次插值;nearest:最鄰近插值法;v4:MATLAB 4中的griddata 算法。命令6 spline功能 三次樣條數(shù)據(jù)插值格式 yy = spline(x,y,xx %對于給定的離散的測量數(shù)據(jù)x,y(稱為斷點(diǎn),要尋找一個(gè)三項(xiàng)多項(xiàng)式x (p y =,以逼近每對數(shù)據(jù)(x,y點(diǎn)間的曲線。過兩點(diǎn)y ,x (i i 和y ,x (1i 1i +只能確定一條直線,而通過一點(diǎn)的三次多項(xiàng)式曲線有無窮多條。為使通過中間斷點(diǎn)的三次多項(xiàng)式曲線具有唯一

52、性,要增加兩個(gè)條件(因?yàn)槿味囗?xiàng)式有4個(gè)系數(shù):1.三次多項(xiàng)式在點(diǎn)y ,x (i i 處有: x (p x (p i i i i ''='2.三次多項(xiàng)式在點(diǎn)y ,x (1i 1i +處有:x (p x (p 1i i 1i i +''=' 3.p(x在點(diǎn)y ,x (i i 處的斜率是連續(xù)的(為了使三次多項(xiàng)式具有良好的解析性,加上的條件;4.p(x在點(diǎn)y ,x (i i 處的曲率是連續(xù)的;對于第一個(gè)和最后一個(gè)多項(xiàng)式,人為地規(guī)定如下條件:.x (p x (p 21'''=''' .x (p x (p 1n

53、n -'''=''' 上述兩個(gè)條件稱為非結(jié)點(diǎn)(not-a-knot條件。綜合上述內(nèi)容,可知對數(shù)據(jù)擬合的三次樣條函數(shù)p(x是一個(gè)分段的三次多項(xiàng)式:=+1n n n 322211x x x x (p x x x x (p x x x x (p x (p ,其中每段x (p i 都是三次多項(xiàng)式。 該命令用三次樣條插值計(jì)算出由向量x 與y 確定的一元函數(shù)y=f(x在點(diǎn)xx 處的值。若參量y 是一矩陣,則以y 的每一列和x 配對,再分別計(jì)算由它們確定的函數(shù)在點(diǎn)xx 處的值。則yy 是一階數(shù)為length(xx*size(y,2的矩陣。pp = spline

54、(x,y %返回由向量x 與y 確定的分段樣條多項(xiàng)式的系數(shù)矩陣pp ,它可用于命令ppval 、unmkpp 的計(jì)算。例2-36對離散地分布在y=exp(xsin(x函數(shù)曲線上的數(shù)據(jù)點(diǎn)進(jìn)行樣條插值計(jì)算:>>x = 0 2 4 5 8 12 12.8 17.2 19.9 20; y = exp(x.*sin(x;>>xx = 0:.25:20;>>yy = spline(x,y,xx;>>plot(x,y,'o',xx,yy 插值圖形結(jié)果為圖2-19。 圖2-19 三次樣條插值命令7 interpn功能n維數(shù)據(jù)插值(查表格式VI =

55、 interpn(X1,X2,Xn,V,Y1,Y2,Yn %返回由參量X1,X2,Xn,V確定的n元函數(shù)V=V(X1,X2,Xn在點(diǎn)(Y1,Y2,Yn處的插值。參量Y1,Y2,Yn是同型的矩陣或向量。若Y1,Y2,Yn是向量,則可以是不同長度,不同方向(行或列的向量。它們將通過命令ndgrid生成同型的矩陣,再作計(jì)算。若點(diǎn)(Y1,Y2,Yn中有位于點(diǎn)(X1,X2,Xn之外的點(diǎn),則相應(yīng)地返回特殊變量NaN。VI = interpn(V,Y1,Y2,Yn %缺省地,X1=1:size(V,1,X2=1:size(V,2,Xn=1:size(V,n,再按上面的情形計(jì)算。VI = interpn(V,

56、ntimes %作ntimes次遞歸計(jì)算,在V的每兩個(gè)元素之間插入它們的n維插值。這樣,V的階數(shù)將不斷增加。interpn(V等價(jià)于interpn(V,1。VI = interpn(,method %用指定的算法method計(jì)算:linear:線性插值(缺省算法;cubic:三次插值;spline:三次樣條插值法;nearest:最鄰近插值算法。命令8 meshgrid功能生成用于畫三維圖形的矩陣數(shù)據(jù)。格式X,Y = meshgrid(x,y 將由向量x,y(可以是不同方向的指定的區(qū)域min(x, max(x,min(y,max(y用直線x=x(i,y=y(j(i=1,2,lengt h(x ,j=1,2,length(y進(jìn)行劃分。這樣,得到了le