Matlab學(xué)習(xí)與實(shí)驗(yàn)教程-第7章 MATLAB數(shù)值積分與數(shù)值微分

第7章MATLAB數(shù)值積分與數(shù)值微分MATLAB數(shù)值積分

MATLAB數(shù)值微分MATLAB離散傅里葉變換7.1數(shù)值積分7.1.1數(shù)值積分基本原理求解定積分的數(shù)值方法多種多樣，如簡(jiǎn)單的梯形法、辛普生(Simpson)法、牛頓－柯特斯(Newton-Cotes)法等都是經(jīng)常采用的方法。它們的基本思想都是將整個(gè)積分區(qū)間[a,b]分成n個(gè)子區(qū)間[xi,xi+1]，i=1,2,…,n，其中x1=a，xn+1=b。這樣求定積分問題就分解為求和問題。7.1.2數(shù)值積分的實(shí)現(xiàn)方法1．自適應(yīng)積分法MATLAB提供了基于全局自適應(yīng)積分算法的integral函數(shù)來求定積分。函數(shù)的調(diào)用格式為：I=integral(@fname,a,b)其中，I是計(jì)算得到的積分；filename是被積函數(shù)；a和b分別是定積分的下限和上限，積分限可以為無窮大。例7-1求。2．變步長辛普生法基于變步長辛普生法，MATLAB給出了quad函數(shù)和quadl函數(shù)來求定積分。函數(shù)的調(diào)用格式為：[I,n]=quad(@fname,a,b,tol,trace)[I,n]=quadl(@fname,a,b,tol,trace)其中，fname是被積函數(shù)名。a和b分別是定積分的下限和上限。tol用來控制積分精度，默認(rèn)時(shí)取tol=10-6。trace控制是否展現(xiàn)積分過程，若取非0則展現(xiàn)積分過程，取0則不展現(xiàn)，默認(rèn)時(shí)取trace=0。返回參數(shù)I即定積分值，n為被積函數(shù)的調(diào)用次數(shù)。

例7-2設(shè)f(x)=，求S=。(1)建立被積函數(shù)文件fesin.m。functionf=fesin(x)f=exp(-0.5*x).*sin(x+pi/6);(2)調(diào)用數(shù)值積分函數(shù)quad求定積分。[S,n]=quad('fesin',0,3*pi)S=0.9008n=77在建立被積函數(shù)文件時(shí)，注意到所寫程序應(yīng)允許向量作為輸入?yún)?shù)，所以在該函數(shù)文件中采用了.*運(yùn)算符。也可不建立關(guān)于被積函數(shù)的函數(shù)文件，而使用語句函數(shù)（內(nèi)聯(lián)函數(shù)）求解，命令如下：g=inline('exp(-0.5*x).*sin(x+pi/6)');%定義一個(gè)語句函數(shù)[S,n]=quad(g,0,3*pi) %注意函數(shù)名不加@號(hào)S=0.9008n=77

例7-3分別用quad函數(shù)和quadl函數(shù)求的近似值，并在相同的積分精度下，比較函數(shù)的調(diào)用次數(shù)。調(diào)用函數(shù)quad求定積分：formatlongfx=inline('exp(-x)');[I,n]=quad(fx,1,2.5,1e-10)I=0.285794442547663n=65調(diào)用函數(shù)quadl求定積分：formatlongfx=inline('exp(-x)');[I,n]=quadl(fx,1,2.5,1e-10)I=

0.285794442548811n=18formatshort3．高斯-克朗羅德法MATLAB提供了基于自適應(yīng)高斯-克朗羅德法的quadgk函數(shù)來求振蕩函數(shù)的定積分。該函數(shù)的調(diào)用格式為[I,err]=quadgk(@fname,a,b)其中，err返回近似誤差范圍，其他參數(shù)的含義和用法與quad函數(shù)相同。積分上下限可以是?Inf或Inf，也可以是復(fù)數(shù)。如果積分上下限是復(fù)數(shù)，則quadgk在復(fù)平面上求積分。

例7-4求定積分。(1)被積函數(shù)文件fx.m。functionf=fx(x)f=x.*sin(x)./(1+cos(x).*cos(x));(2)調(diào)用函數(shù)quad8求定積分。I=quad8('fx',0,pi)I=2.46744．梯形積分法在科學(xué)實(shí)驗(yàn)和工程應(yīng)用中，函數(shù)關(guān)系往往是不知道的，只有實(shí)驗(yàn)測(cè)定的一組樣本點(diǎn)和樣本值，這時(shí)，人們就無法使用quad等函數(shù)計(jì)算其定積分。在MATLAB中，對(duì)由表格形式定義的函數(shù)關(guān)系的求定積分問題用梯形積分函數(shù)trapz。該函數(shù)調(diào)用格式如下：I=trapz(X,Y)其中，向量X，Y定義函數(shù)關(guān)系Y

=

f(X)。

例7-5用trapz函數(shù)計(jì)算定積分。命令如下：X=1:0.01:2.5;Y=exp(-X);%生成函數(shù)關(guān)系數(shù)據(jù)向量trapz(X,Y)ans=0.2858

7.1.3多重定積分的數(shù)值求解MATLAB提供的dblquad函數(shù)用于求二重積分的數(shù)值解，triplequad函數(shù)用于求三重積分的數(shù)值解。函數(shù)的調(diào)用格式為dblquad(@fun,a,b,c,d,tol)triplequad(@fun,a,b,c,d,e,f,tol)其中，fun為被積函數(shù)，[a，b]為x的積分區(qū)域，[c，d]為y的積分區(qū)域，[e，f

]為z的積分區(qū)域，參數(shù)tol的用法與函數(shù)quad完全相同。

例7-5計(jì)算二重定積分

（1）建立一個(gè)函數(shù)文件fxy.m：

functionf=fxy(x,y)

globalki;

ki=ki+1;%ki用于統(tǒng)計(jì)被積函數(shù)的調(diào)用次數(shù)

f=exp(-x.^2/2).*sin(x.^2+y);

（2）調(diào)用dblquad函數(shù)求解。

globalki;

ki=0;

I=dblquad(@fxy,-2,2,-1,1)

ki

I=

1.5745

ki=

1050如果使用匿名函數(shù)，則命令如下：

f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');

I=dblquad(f,-2,2,-1,1)

I=

1.5745例7-6計(jì)算三重定積分

命令如下：

fxyz=inline('4*x.*z.*exp(-z.*z.*y-x.*x)','x','y','z');

triplequad(fxyz,0,pi,0,pi,0,1,1e-7)

ans=

1.73287.2數(shù)值微分7.2.1數(shù)值差分與差商7.2.2數(shù)值微分的實(shí)現(xiàn)在MATLAB中，沒有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，只有計(jì)算向前差分的函數(shù)diff，其調(diào)用格式為：DX=diff(X)：計(jì)算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,…,n-1。DX=diff(X,n)：計(jì)算X的n階向前差分。例如，diff(X,2)=diff(diff(X))。DX=diff(A,n,dim)：計(jì)算矩陣A的n階差分，dim=1時(shí)(缺省狀態(tài))，按列計(jì)算差分；dim=2，按行計(jì)算差分。例7-6生成以向量V=[1,2,3,4,5,6]為基礎(chǔ)的范得蒙矩陣，按列進(jìn)行差分運(yùn)算。命令如下：V=vander(1:6)DV=diff(V)%計(jì)算V的一階差分例7-7用不同的方法求函數(shù)f(x)的數(shù)值導(dǎo)數(shù)，并在同一個(gè)坐標(biāo)系中做出f'(x)的圖像。程序如下：f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2');g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5');x=-3:0.01:3;p=polyfit(x,f(x),5);%用5次多項(xiàng)式p擬合f(x)dp=polyder(p);%對(duì)擬合多項(xiàng)式p求導(dǎo)數(shù)dpdpx=polyval(dp,x);%求dp在假設(shè)點(diǎn)的函數(shù)值dx=diff(f([x,3.01]))/0.01;%直接對(duì)f(x)求數(shù)值導(dǎo)數(shù)gx=g(x);%求函數(shù)f的導(dǎo)函數(shù)g在假設(shè)點(diǎn)的導(dǎo)數(shù)plot(x,dpx,x,dx,'.',x,gx,'-');%作圖7.3離散傅里葉變換離散傅里葉變換（DFT）廣泛應(yīng)用于信號(hào)分析、光譜和聲譜分析、全息技術(shù)等各個(gè)領(lǐng)域。但直接計(jì)算DFT的運(yùn)算量與變換的長度N的平方成正比，當(dāng)N較大時(shí)，計(jì)算量太大。隨著計(jì)算機(jī)技術(shù)的迅速發(fā)展，在計(jì)算機(jī)上進(jìn)行離散傅里葉變換計(jì)算成為可能，特別是快速傅里葉變換（FFT）算法的出現(xiàn)，為離散傅里葉變換的應(yīng)用創(chuàng)造了條件。MATLAB提供了一套計(jì)算快速傅里葉變換的函數(shù)，它們包括求一維、二維和N維離散傅里葉變換函數(shù)fft、fft2和fftn，還包括求上述各維離散傅里葉變換的逆變換函數(shù)ifft、ifft2和ifftn等。本節(jié)先簡(jiǎn)要介紹離散傅里葉變換的基本概念和變換公式，然后討論MATLAB中離散傅里葉變換的實(shí)現(xiàn)。7.3.1離散傅里葉變換算法簡(jiǎn)介在某時(shí)間片等距地抽取N個(gè)抽樣時(shí)間tm處的樣本值f(tm)，且記作f(m)，這里m=0，1，2，…，N-1，稱向量F(k)（k=0，1，2，…，N-1）為f(m)的一個(gè)離散傅里葉變換，其中因?yàn)镸ATLAB不允許有零下標(biāo)，所以將上述公式中m的下標(biāo)均移動(dòng)1，于是便得相應(yīng)公式：由f(m)求F(k)的過程，稱為求f(m)的離散傅里葉變換，或稱為F(k)為f(m)的離散頻譜。反之，由F(k)逆求f(m)的過程，稱為離散傅里葉逆變換，相應(yīng)的變換公式為：7.3.2離散傅里葉變換的實(shí)現(xiàn)MATLAB提供了對(duì)向量或直接對(duì)矩陣進(jìn)行離散傅里葉變換的函數(shù)。下面只介紹一維離散傅里葉變換函數(shù)，其調(diào)用格式與功能為：（1）fft(X)：返回向量X的離散傅里葉變換。設(shè)X的長度（即元素個(gè)數(shù)）為N，若N為2的冪次，則為以2為基數(shù)的快速傅里葉變換，否則為運(yùn)算速度很慢的非2冪次的算法。對(duì)于矩陣X，fft(X)應(yīng)用于矩陣的每一列。（2）fft(X,N)：計(jì)算N點(diǎn)離散傅里葉變換。它限定向量的長度為N，若X的長度小于N，則不足部分補(bǔ)上零；若大于N，則刪去超出N的那些元素。對(duì)于矩陣X，它同樣應(yīng)用于矩陣的每一列，只是限定了向量的長度為N。（3）fft(X,[],dim)或fft(X,N,dim)：這是對(duì)于矩陣而言的函數(shù)調(diào)用格式，前者的功能與FFT(X)基本相同，而后者則與FFT(X,N)基本相同。只是當(dāng)參數(shù)dim=1時(shí)，該函數(shù)