第6章MATLAB數(shù)值計算

1、第第6 6章章MATLABMATLAB數(shù)值計算數(shù)值計算6.1 6.1 數(shù)據(jù)處理與多項式計算數(shù)據(jù)處理與多項式計算6.2 6.2 傅立葉分析傅立葉分析6.3 6.3 數(shù)值微積分?jǐn)?shù)值微積分6.4 6.4 線性方程組求解線性方程組求解 6.5 6.5 常微分方程的數(shù)值求解常微分方程的數(shù)值求解6.6 6.6 非線性方程的數(shù)值求解非線性方程的數(shù)值求解6.7 6.7 稀疏矩陣稀疏矩陣6.1 數(shù)據(jù)處理與多項式計算數(shù)據(jù)處理與多項式計算6.1.1 數(shù)據(jù)統(tǒng)計與分析數(shù)據(jù)統(tǒng)計與分析1. 求矩陣最大和最小元素(1)求向量的最大最小元素y=max(X) 返回向量X的最大元素存入y。y,I=max(X) 返回向量X的最大元素

2、存入y，最大元素的序號存入I。(2)求矩陣的最大和最小元素max(A) 返回一個行向量行向量，向量的第i個元素是A矩陣的第i列列上的最大元素。Y,U=max(A) 返回兩個行向量，Y向量記錄A的每列列的最大元素，U向量記錄每列最大元素的行行號。max(A,dim) dim取1或2。dim取1時，該函數(shù)和max(A)完全相同。dim取2時，該函數(shù)返回一個列向量列向量，其第i個元素是A矩陣的第i行上的最大元素。(3)兩個向量或矩陣對應(yīng)元素的比較U=max(A,B) A,B是兩個同型的向量或矩陣。結(jié)果U是與A,B同型的向量或矩陣，U的每個元素等于A,B對應(yīng)元素的較大者。U=max(A,n) n是一個

3、標(biāo)量。結(jié)果U是與A同型的向量或矩陣，U的每個元素等于A對應(yīng)元素和n中的較大者。min函數(shù)的用法和max完全相同。例例6.1 求矩陣A的每行及每列的最大和最小元素，并求整個矩陣的最大和最小元素。命令如下：A=13,-56,78;25,63,-235;78,25,563;1,0,-1;max(A,2) %求每行最大元素求每行最大元素min(A,2) %求每行最小元素求每行最小元素max(A) %求每列最大元素求每列最大元素min(A) %求每列最小元素求每列最小元素 %求整個矩陣的最大元素求整個矩陣的最大元素 %求整個矩陣的最小元素求整個矩陣的最小元素max(max(A)min(min(A)2.

4、求矩陣的平均值和中值 求矩陣和向量元素的平均值的函數(shù)是mean ; 求中值的函數(shù)是median。 它們的調(diào)用方法和max函數(shù)完全相同。3. 矩陣元素求和與求積 矩陣和向量求和與求積的基本函數(shù)是sum和prod， 其使用方法和max類似。例例6.2 求矩陣A的每行元素的乘積和全部元素的乘積。命令如下：A=1,2,3,4;5,6,7,8;9,10,11,12;S=prod(A,2) %每行元素乘積prod(S) %求A的全部元素的乘積（每列元素的積）4. 矩陣元素累加和與累乘積MATLAB中，使用cumsum和cumprod函數(shù)能方便地求得向量和矩陣元素的累加和與累乘積向量，函數(shù)的用法和sum及p

5、rod相同例例6.3求向量X=(1！，2！，3！，10！)。命令如下：X=cumprod(1:10)5. 標(biāo)準(zhǔn)方差MATLAB中，提供了計算數(shù)據(jù)序列的標(biāo)準(zhǔn)方差的函數(shù)std。對于向量X，std(X)返回一個標(biāo)準(zhǔn)方差。對于矩陣A，std(A)返回一個行向量，它的各個元素便是矩陣A各列或各行的標(biāo)準(zhǔn)方差。 std函數(shù)的一般調(diào)用格式為：std(A,FLAG,dim) 其中dim取1或2。 當(dāng)dim=1時，求各列元素的標(biāo)準(zhǔn)方差； 當(dāng)dim=2時，則求各行元素的標(biāo)準(zhǔn)方差。 FLAG取0或1。6. 元素排序MATLAB中對向量X是排序函數(shù)是sort(X)，函數(shù)返回一個對X中的元素按升序升序排列的新向量。sor

6、t函數(shù)也可以對矩陣A的各列(或行)重新排序，其調(diào)用格式為：Y,I=sort(A,dim) 其中dim指明對A的列還是行進行排序，若dim=1，則按列排，若dim=2，則按行排。Y是排序后的矩陣，而I記錄Y中的元素在A中位置。例例6.4 對矩陣做各種排序。命令如下：A=1,-8,5;4,12,6;13,7,-13;sort(A) %對A的每列按升序排序-sort(-A,2) %對A的每行按降序排序X,I=sort(A) %對A按列排序，并將每個元素所在 行號送矩陣I6.1.2 數(shù)值插值數(shù)值插值1. 一維數(shù)值插值interp1函數(shù)調(diào)用格式為：Y1=interp1(X,Y,X1,method)函數(shù)根

7、據(jù)X、Y的值，計算函數(shù)在X1處的值。X、Y是兩個等長的已知向量，分別描述采樣點和樣本值，X1是一個向量或標(biāo)量，描述欲插值的點，Y1是一個與X1等長的插值結(jié)果。method是插值方法，允許的取值有l(wèi)inear(線性插值)、nearest(最近插值)、spline(三次樣條插值)、cubic（三次多項式插值），缺省值是linear。例例6.5 用不同的插值方法計算sin(x)在/2點的值。這是一個一維插值問題。在MATLAB命令窗口，輸入命令：X=0:0.2:pi;Y=sin(X); %給出X、Yinterp1(X,Y,pi/2) %用缺省方法(即線性插值方法)計算 sin(/2)interp1(

8、X,Y,pi/2,nearest) %用最近方法計算sin(/2)interp1(X,Y,pi/2,linear) %用線性方法計算sin(/2)interp1(X,Y,pi/2,spline) %用三次樣條方法計算sin(/2)interp1(X,Y,pi/2,cubic) %用三次多項式方法計算sin(/2)M A T L A B 中 有 一 個 專 門 的 三 次 樣 條 插 值 函 數(shù)Y1=spline(X,Y,X1)，其功能及使用方法與函數(shù)Y1=interp1(X,Y,X1,spline)完全相同。 例例6.6 已知檢測參數(shù)f隨時間t的采樣結(jié)果，用數(shù)值插值法計算t=2，7，12，17

9、，22，17，32，37，42，47，52，57時f的值。這是一個一維數(shù)值插值問題，命令如下：T=0:5:65;X=2:5:57;F=3.2015,2.2560,879.5,1835.9,2968.8,4136.2,5237.9,6152.7,.6725.3,6848.3,6403.5,6824.7,7328.5,7857.6;F1=interp1(T,F,X) %用線性方法插值F1=interp1(T,F,X,nearest) %用最近方法插值F1=interp1(T,F,X,spline) %用三次樣條方法插值F1=interp1(T,F,X,cubic) %用三次多項式方法插值2. 二維

10、數(shù)值插值MATLAB中，提供了解決二維插值問題的函數(shù)。其調(diào)用格式為：Z1=interp2(X,Y,Z,X1,Y1,method) 其中，X、Y是兩個向量，分別描述兩個參數(shù)的采樣點，Z是與參數(shù)采樣點對應(yīng)的采樣變量的樣本值，X1、Y1是兩個向量或標(biāo)量，描述欲插值的點。method的取值與一維插值函數(shù)相同。例例6.7a 設(shè)Z=x2+y2，對Z函數(shù)在(0，1)(0，2)區(qū)域內(nèi)進行插值。命令如下：x=0:0.1:10;y=0:0.2:20;X,Y=meshgrid(x,y);Z=X.2+Y.2;interp2(x,y,Z,0.5,0.5) %對函數(shù)在(0.5,0.5)點進行插值interp2(x,y,Z

11、,0.5 0.6,0.4) %對函數(shù)在(0.5,0.4)點和(0.6，0.4) 點進行插值interp2(x,y,Z,0.5 0.6,0.4 0.5) %對函數(shù)在(0.5,0.4)點和(0.6， 0.5)點進行插值interp2(x,y,Z,0.5 0.6 ,0.4 0.5) %對函數(shù)在 (0.5,0.4),(0.6,0.4),(0.5,0.5)和 (0.6,0.5)點進行插值例6.7b 某實驗對一根長10米的鋼軌進行熱源的溫度傳播測試。用x表示測量點0:2.5:10(米)，用h表示測量時間0:30:60(秒)，用T表示測試所得各點的溫度()。試用線性插值求出在一分鐘內(nèi)每隔10秒、鋼軌每隔0.

12、5米處的溫度。程序如下：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:0.5:10;hi=0:10:60;temps=interp2(x,h,T,xi,hi,cubic);mesh(xi,hi,temps);02468100204060-200204060801003. 三維數(shù)值插值對三維函數(shù)插值的函數(shù)是interp3，其使用方法和interp2相同。其調(diào)用格式為：W1=interp3(X,Y,Z,W,X1,Y1,Z1,method)函數(shù)返回三維插值結(jié)果。 其中，X、Y、Z是三個向量，分別描述三個參數(shù)的

13、采樣點，W是與參數(shù)采樣點對應(yīng)的采樣變量的樣本值，X1、Y1、Z1是三個向量或標(biāo)量，描述欲插值的點。method是插值方法，可選，其缺省值是 line。method的取值與一、二維插值函數(shù)相同。6.1.3 曲線擬合曲線擬合MATLAB中，提供了解決使用最小二乘法進行曲線擬合的函數(shù)。調(diào)用格式為：P,S=polyfit(X,Y,m) 函數(shù)根據(jù)采樣點X和采樣點函數(shù)值Y，產(chǎn)生一個m次多項式P及其在采樣點的誤差向量S。 其中X、Y是兩個等長的向量，P是一個長度為m+1的向量。例例6.8 用一個3次多項式在區(qū)間0，2內(nèi)逼近函數(shù)sin(x)。命令如下：X=linspace(0,2*pi,50);Y=sin(X

14、);P,S=polyfit(X,Y,3)plot(X,Y, -*,X,polyval(P,X), :o)01234567-1-0.8-0.6-0.4-0.200.20.40.60.816.1.4 多項式計算多項式計算1. 多項式的建立已知一個多項式的全部根X求多項式系數(shù)的函數(shù)是poly(X)，該函數(shù)返回以X為全部根的一個多項式P，當(dāng)X是一個長度為m的向量時，P是一個長度為m+1的向量。2. 多項式求根求多項式p(x)的根的函數(shù)是roots(P)，這里，P是p(x)的系數(shù)向量，該函數(shù)返回方程p(x)=0的全部根(含重根，復(fù)根)。3. 多項式求值求多項式p(x)在某點或某些點的函數(shù)值的函數(shù)是pol

15、yval(P,x)。若x為一數(shù)值，則求多項式在該點的值；若x為向量或矩陣，則對向量或矩陣中的每個元素求其多項式的值。例例6.9 已知一個多項式，計算： (1)計算f(x)=0 的全部根。(2)由方程f(x)=0的根構(gòu)造一個多項式g(x)，并與f(x)進行對比。(3)計算f(5)、f(7.8)、f(9.6)、f(12.3)的值。命令如下：P=3,0,4,-5,-7.2,5;X=roots(P) %求方程f(x)=0的根G=poly(X) %求多項式g(x)X0=5,7.8,9.6,12.3;f=polyval(P,X0) %求多項式f(x)在給定點的值多項式求值還有一個函數(shù)是polyvalm，其

16、調(diào)用格式與polyval相同，但含義不同。polyvalm函數(shù)要求x為方陣，它以方陣為自變量求多項式的值。 4. 多項式的四則運算(1)多項式的加減法當(dāng)兩個多項式的次數(shù)不同時，要在一個較低次冪的多項式系數(shù)向量前補0，使兩個系數(shù)向量等長。(2)多項式的乘法函數(shù)conv(P1,P2)用于求多項式P1和P2的乘積。(3)多項式的除法函數(shù)Q,r=deconv(P1,P2)用于對多項式P1和P2作除法運算。其中Q返回多項式P1除以P2的商式，r返回P1除以P2的余式。這里，Q和r仍是多項式系數(shù)向量。deconv是conv的逆函數(shù)，即有P1=conv(P2,Q)+r。例例6.10 設(shè)有兩個多項式，計算：(

17、1)求f(x)+g(x)、f(x)-g(x)。(2)求f(x)g(x)、f(x)/g(x)。在MATLAB命令窗口，輸入命令：f=3,-5,2,-7,5,6;g=3,5,-3;g1=0,0,0,g;f+g1 %求f(x)+g(x)f-g1 %求f(x)-g(x)conv(f,g) %求f(x)*g(x)Q,r=deconv(f,g) %求f(x)/g(x)，商式送Q，余式送r。5. 多項式的導(dǎo)函數(shù)對多項式求導(dǎo)數(shù)的函數(shù)是：p=polyder(P) 求多項式P的導(dǎo)函數(shù)p=polyder(P,Q) 求P*Q的導(dǎo)函數(shù)p,q=polyder(P,Q) 求P/Q的導(dǎo)函數(shù)，導(dǎo)函數(shù)的分子存入p，分母存入q。例

18、例6.11 求有理分式的導(dǎo)數(shù)。命令如下：P=3,5,0,-8,1,-5;Q=10,5,0,0,6,0,0,7,-1,0,-100;p,q=polyder(P,Q)6.1.5 函數(shù)的最大值與最小值函數(shù)的最大值與最小值MATLAB中用于求最小值的函數(shù)是：fmin(f,a,b) 求單變量函數(shù)f(x)在區(qū)間(a,b)上的最小值點。fmins(F,X0) 求多變量函數(shù)F(x)在估計值X0附近的最小值點。MATLAB沒有專門提供求函數(shù)最大值點的函數(shù)，但只要注意到-f(x)在區(qū)間(a,b)上的最小值點就是f(x)在(a,b)的最大值點，所以fmin(-f,a,b)返回函數(shù)f(x)在區(qū)間(a,b)上的最大值。

19、例例6.12 求函數(shù)f(x)在區(qū)間(-10，1)和(1，10)上的最小值點。首先建立函數(shù)文件fx.m：function f=f(x)f=x-1/x+5;return上述函數(shù)文件也可用一個語句函數(shù)代替：Ff=inline( x-1/x+5 )再在MATLAB命令窗口，輸入命令：fmin(fx,-10,-1) %求函數(shù)在區(qū)間(-10，-1)內(nèi)的最小值點fmin(ff,1,10) %求函數(shù)在區(qū)間(1，10)內(nèi)的最小值點。 注意函數(shù)名f不用加有時為了描述某個數(shù)學(xué)函數(shù)的方便，可以用inline()函數(shù)來直接編寫該函數(shù)，形式相當(dāng)于M-函數(shù)，但無編寫一個真正的MATLAB文件，就可以描述出某種數(shù)學(xué)關(guān)系。其調(diào)

20、用格式為fun=inline(函數(shù)內(nèi)容,自變量列表)例例6.13 設(shè)有函數(shù)f(x,y,z)，求函數(shù)f在(0.5,0.5,0.5)附近的最小值。建立函數(shù)文件fxyz.m：function f=f(u)x=u(1);y=u(2);z=u(3);f=x+y.2./x/4+z.2./y+2./z;return在MALAB命令窗口，輸入命令：U=fmins(fxyz,0.5,0.5,0.5) %求函數(shù)的最小值點fxyz(U) %求函數(shù)的最小值6.2 6.2 傅立葉分析傅立葉分析MATLAB中，提供了對向量(或直接對矩陣的行或列)進行離散傅立葉變換的函數(shù)，其調(diào)用格式是：Y=fft(X,n,dim)(1)當(dāng)

21、X是一個向量時，返回對X的離散傅立葉變換。(2)當(dāng)X是一個矩陣時，返回一個矩陣并送Y，其列(行)是對X的列(行)的離散傅立葉變換。例例6.14 求X=(1，0，-3，5，2)的離散傅立葉變換。在MATLAB命令窗口，輸入命令：X=1,0,-3,5,2;Y=fft(X) %對X進行變換3. 離散傅立葉變換的逆變換MATLAB中，對向量(或直接對矩陣的行或列)進行離散傅立葉逆變換的函數(shù)的調(diào)用方法是：Y=ifft(X,n,dim)函數(shù)對X進行離散傅立葉逆變換。其中X、n、dim的意義及用法和離散傅立葉變換函數(shù)fft完全相同。例例6.15 對矩陣A的列向量、行向量分別進行離散傅立葉變換、并對變換結(jié)果進

22、行逆變換。命令如下：A=3,2,1,1;-5,1,0,1;3,2,1,5;fftA=fft(A) %求A的列向量的傅立葉變換fftA2=fft(A,4,2) %求A的行向量的傅立葉變換ifft(fftA) %對矩陣fftA的列向量進行傅立葉逆 變換，結(jié)果應(yīng)等于Aifft(fftA2,4,2) %對矩陣fftA2的行向量進行傅立 葉逆變換，其結(jié)果應(yīng)等于AfftA = 1.0000 5.0000 2.0000 7.0000 4.0000 + 6.9282i 0.5000 + 0.8660i 0.5000 + 0.8660i -2.0000 + 3.4641i 4.0000 - 6.9282i 0.

23、5000 - 0.8660i 0.5000 - 0.8660i -2.0000 - 3.4641ifftA2 = 7.0000 2.0000 - 1.0000i 1.0000 2.0000 + 1.0000i -3.0000 -5.0000 -7.0000 -5.0000 11.0000 2.0000 + 3.0000i -3.0000 2.0000 - 3.0000ians = 3.0000 2.0000 1.0000 1.0000 -5.0000 1.0000 0.0000 1.0000 3.0000 2.0000 1.0000 5.0000ans = 3 2 1 1 -5 1 0 1 3

24、 2 1 56.3 6.3 數(shù)值微積分?jǐn)?shù)值微積分6.3.1 數(shù)值微分MATLAB中，沒有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，只有計算向前差分的函數(shù)。DX=diff(X) 計算向量X的向前差分， DX(i)=X(i+1)-X(i)，0i0 %非齊次方程組 if rank(A)=rank(a,b) %方程組相容 if rank(A)=m %有唯一解 x=Ab; else %方程組有無窮多個解,基礎(chǔ)解系 disp(原方程組有有無窮個解，其齊次方程組的基礎(chǔ)解系為y，特解為x); y=null(A,r); x=Ab; end else %方程組不相容，給出最小二乘法解 disp(方程組的最小二乘法解是：); x=

25、Ab; end else %齊次方程組 if rank(A)=n %列滿秩 x=zero(m,1) %0解 else %非0解 disp(方程組有無窮個解，基礎(chǔ)解系為x); x=null(A,r); end endreturn2. 應(yīng)用舉例例例6.24 求線性方程組的解。在MATLAB命令窗口，輸入命令：A=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;b=4,6,12,6;x,y=line_solution(A,b) %調(diào)用自定義函數(shù)例例6.25 求下列線性方程組的解。在MATLAB命令窗口，輸入命令：A=2,7,3,1;3,5,2,2;9,4,1,7;b=6,4,

26、2;x,y=line_solution(A,b)6.5 6.5 常微分方程的數(shù)值求解常微分方程的數(shù)值求解基于龍格庫塔法，MATLAB提供了求常微分方程數(shù)值解的函數(shù)，一般調(diào)用格式為： X,Y=ode23(f,x0,xn,y0) X,Y=ode45(f,x0,xn,y0) 其中，X、Y是兩個向量， X對應(yīng)自變量x在求解區(qū)間x1，xn的一組采樣點，其采樣密度是自適應(yīng)的，無需指定； Y是與X對應(yīng)的一組解， f是一個函數(shù)， x0，xn代表自變量的求解區(qū)間，y0=y(x0)，由方程的初值給定。 函數(shù)在求解區(qū)間x0，xn內(nèi)，自動設(shè)立自動設(shè)立采樣點向量X，并求出解函數(shù)y在采樣點X處的樣本值。 例例6.26 求

27、微分方程初值問題在1，3區(qū)間內(nèi)的數(shù)值解，并將結(jié)果與解析解進行比較。先建立一個該函數(shù)的m文件fxy1.m：function f=f(x,y)f=-2*y./x+4*x %注意使用點運算符return再輸入命令：X,Y=ode45(fxy1,1,3,2);X %顯示自變量的一組采樣點Y %顯示求解函數(shù)與采樣點對應(yīng)的一組數(shù)值解(X.2+1./X.2) %顯示求解函數(shù)與采樣點對應(yīng)的一組解析解例例6.27 求解初值問題在區(qū)間0，2中的解。建立一個函數(shù)文件 fxy2.m：function f=f(x,y)f(2)=-x.*y(2)+x.2-5;f(1)=y(2);f=f;return在MATLAB命令窗口

28、，輸入命令：X,Y=ode45(fxy2,0,2,5,6);X,Y6.6 6.6 非線性方程的數(shù)值求解非線性方程的數(shù)值求解1單變量非線性方程求解M AT L A B 中 ， 提 供 了 求 解 單 變 量 方 程 的 函 數(shù)fzero(f,x0,tol)，該函數(shù)采用迭代法計算函數(shù)f(x)的一個零點，迭代初值為x0，當(dāng)兩次迭代結(jié)果小于tol時停止迭代過程。tol的缺省值是eps。 注意，在調(diào)用函數(shù)fzero 之前，要使用m文件建立自己要計算的函數(shù)f(x)，只有定義了函數(shù)f(x)的m文件后，才能在fzero函數(shù)的參數(shù)中使用自定義函數(shù)名。例例6.28 求f(x)=x-1/x+5 在x0=-5和x0=

29、1作為迭代初值時的零點。先編制一個函數(shù)文件fz.m：function f=f(x)f=x-1/x+5;然后，在MATLAB命令窗口，輸入命令：fzero(fz,-5) %以-5作為迭代初值Zero found in the interval: -4.8, -5.2.fzero(fz,1)2. 非線性方程組的求解非線性方程組的求解 對于非線性方程組對于非線性方程組F(X)=0，用，用fsolve函數(shù)求其數(shù)值函數(shù)求其數(shù)值解。解。fsolve函數(shù)的調(diào)用格式為：函數(shù)的調(diào)用格式為： X=fsolve(fun,X0,option) 其中其中X為返回的解，為返回的解，fun是用于定義需求解的非線性是用于定義

30、需求解的非線性方程組的函數(shù)文件名，方程組的函數(shù)文件名，X0是求根過程的初值，是求根過程的初值，option為最優(yōu)化工具箱的選項設(shè)定。最優(yōu)化工具箱提供了為最優(yōu)化工具箱的選項設(shè)定。最優(yōu)化工具箱提供了20多個選項，用戶可以使用多個選項，用戶可以使用optimset命令將它們顯示出命令將它們顯示出來。如果想改變其中某個選項，則可以調(diào)用來。如果想改變其中某個選項，則可以調(diào)用optimset()函數(shù)來完成。例如，函數(shù)來完成。例如，Display選項決定函數(shù)調(diào)用時中間選項決定函數(shù)調(diào)用時中間結(jié)果的顯示方式，其中結(jié)果的顯示方式，其中off為不顯示，為不顯示，iter表示表示每 步 都 顯 示 ，每 步 都 顯

31、示 ， f i n a l 只 顯 示 最 終 結(jié) 果 。只 顯 示 最 終 結(jié) 果 。optimset(Display,off)將設(shè)定將設(shè)定Display選項為選項為off。例例6.29 求方程組在(1，1，1)附近的解并對結(jié)果進行驗證。首先建立方程的函數(shù)文件fxyz1.m：function F=F(X)x=X(1);y=X(2);z=X(3);F(1)=sin(x)+y+z2*exp(x);F(2)=x+y*z;F(3)=x*y*z;在MATLAB命令窗口，輸入命令：X=fsolve(fxyz1,1,1,1) %求解X的三個分量x、y、zY=fxyz1(X) %檢驗所求結(jié)果X是否滿足原方程組norm(Y) %求Y向量的模例例6.30 求圓和直線的兩個交點。建立方程組函數(shù)文件fxyz2.m：function F=F(X)x=X(1);y=X(2);z=X