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

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

7.1數(shù)值微分7.2數(shù)值積分7.3離散傅里葉變換第7章MATLAB數(shù)值微分與積分7.1數(shù)值微分7.1.1數(shù)值差分與差商任意函數(shù)f(x)在x點的導(dǎo)數(shù)是通過極限定義的：如果去掉上述等式右端的h→0的極限過程，并引進記號：稱△f(x)、▽f(x)及δf(x)分別為函數(shù)在x點處以h（h>0）為步長的向前差分、向后差分和中心差分。第7章MATLAB數(shù)值微分與積分當(dāng)步長h充分小時，稱△f(x)/h、▽f(x)/h及δf(x)/h分別為函數(shù)在x點處以h（h>0）為步長的向前差商、向后差商和中心差商。函數(shù)f在點x的微分接近于函數(shù)在該點的任意種差分，而f在點x的導(dǎo)數(shù)接近于函數(shù)在該點的任意種差商。第7章MATLAB數(shù)值微分與積分7.1.2數(shù)值微分的實現(xiàn)有兩種方式計算任意函數(shù)f(x)在給定點x的數(shù)值導(dǎo)數(shù)：第一種方式是用多項式或樣條函數(shù)g(x)對f(x)進行逼近（插值或擬合），然后用逼近函數(shù)g(x)在點x處的導(dǎo)數(shù)作為f(x)在點x處的導(dǎo)數(shù)；第二種方式是用f(x)在點x處的某種差商作為其導(dǎo)數(shù)。第7章MATLAB數(shù)值微分與積分在MATLAB中，沒有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，只有計算向前差分的函數(shù)diff，其調(diào)用格式為：①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時（默認(rèn)狀態(tài)），按列計算差分；dim=2，按行計算差分。第7章MATLAB數(shù)值微分與積分例7-1設(shè)x由[0，2π]間均勻分布的6個點組成，求sinx的1~3階差分。命令如下：>>X=linspace(0,2*pi,6);>>Y=sin(X)Y=00.95110.5878-0.5878-0.9511-0.0000>>DY=diff(Y)%計算Y的一階差分DY=0.9511-0.3633-1.1756-0.36330.9511>>D2Y=diff(Y,2)%計算Y的二階差分，也可用命令diff(DY)計算D2Y=-1.3143-0.81230.81231.3143>>D3Y=diff(Y,3)%計算Y的三階差分，也可用diff(D2Y)或diff(DY,2)D3Y=0.50201.62460.5020第7章MATLAB數(shù)值微分與積分例7-2設(shè)用不同的方法求函數(shù)f(x)的數(shù)值導(dǎo)數(shù)，并在同一個坐標(biāo)系中做出f‘(x)的圖像。f=@(x)sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2;g=@(x)(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次多項式p擬合f(x)dp=polyder(p); %對擬合多項式p求導(dǎo)數(shù)dpdpx=polyval(dp,x); %求dp在假設(shè)點的函數(shù)值dx=diff(f([x,3.01]))/0.01; %直接對f(x)求數(shù)值導(dǎo)數(shù)gx=g(x); %求函數(shù)f的導(dǎo)函數(shù)g在假設(shè)點的導(dǎo)數(shù)plot(x,dpx,x,dx,'.',x,gx,'-') %作圖第7章MATLAB數(shù)值微分與積分結(jié)果表明，用3種方法求得的數(shù)值導(dǎo)數(shù)比較接近。第7章MATLAB數(shù)值微分與積分7.2數(shù)值積分7.2.1數(shù)值積分基本原理設(shè)：

從高等數(shù)學(xué)中知道，當(dāng)|f(x)-p(x)|<ε時，|I1-I2|<ε(b-a)。這說明，當(dāng)ε充分小時，可用I2近似地代替I1。所以，求任意函數(shù)f(x)在[a，b]上的定積分時，要是難以使用解析的方法求出f(x)的原函數(shù)，則可以尋找一個在[a，b]上與f(x)逼近，但形式上卻簡單且易于求積分的函數(shù)p(x)，用p(x)在[a，b]上的積分值近似地代替f(x)在[a，b]上的積分值。一般選擇被積函數(shù)的插值多項式充當(dāng)這樣的替代函數(shù)。選擇的插值多項式的次數(shù)不同，就形成了不同的數(shù)值積分公式。選擇一次多項式時，稱為梯形公式，選擇二次多項式時，稱為辛普森（Simpson）公式。第7章MATLAB數(shù)值微分與積分基本梯形公式：如果把積分區(qū)間[a，b]分劃為n個等長的子區(qū)間：在每個子區(qū)間[ai，ai+1]上用f(x)的插值多項式p(x)代替f(x)，其逼近效果一般將會比在整個區(qū)間上使用一個統(tǒng)一的插值多項式時更好。這樣就形成了數(shù)值積分復(fù)合公式。對被積函數(shù)f(x)采用一、二次多項式插值，然后對插值多項式求積分，就得到了幾個常見的數(shù)值積分公式：基本辛普森公式：復(fù)合梯形公式：復(fù)合辛普森公式：第7章MATLAB數(shù)值微分與積分7.2.2數(shù)值積分的實現(xiàn)基于變步長辛普森法，MATLAB給出了quad函數(shù)和quadl函數(shù)來求定積分。函數(shù)的調(diào)用格式為：[I,n]=quad(filename,a,b,tol,trace)[I,n]=quadl(filename,a,b,tol,trace)其中，filename是被積函數(shù)名；a和b分別是定積分的下限和上限，積分限[a，b]必須是有限的，不能為無窮大（Inf）；tol用來控制積分精度，默認(rèn)時取tol=10-6；trace控制是否展現(xiàn)積分過程，若取非0則展現(xiàn)積分過程，取0則不展現(xiàn)，默認(rèn)時取trace=0；返回參數(shù)I即定積分值，n為被積函數(shù)的調(diào)用次數(shù)。第7章MATLAB數(shù)值微分與積分例7-3求。先建立一個函數(shù)文件fex.m。functionf=fex(x)f=exp(-x.^2);在建立被積函數(shù)的函數(shù)文件時，因為函數(shù)文件允許向量作為輸入?yún)?shù)，所以在函數(shù)表達式中要采用點運算符。第7章MATLAB數(shù)值微分與積分接下來調(diào)用數(shù)值積分函數(shù)quad求定積分，命令如下：>>[I,n]=quad(@fex,0,1)I=0.7468n=13也可不建立關(guān)于被積函數(shù)的函數(shù)文件，而使用匿名函數(shù)（或內(nèi)聯(lián)函數(shù)）求解，命令如下：>>f=@(x)exp(-x.^2); %用匿名函數(shù)f(x)定義被積函數(shù)>>[I,n]=quad(f,0,1)%注意函數(shù)名不加@號I=0.7468n=13第7章MATLAB數(shù)值微分與積分例7-4分別用quad函數(shù)和quadl函數(shù)求的近似值，并在相同的積分精度下，比較函數(shù)的調(diào)用次數(shù)。>>formatlong>>f=@(x)4./(1+x.^2);>>[I,n]=quad(f,0,1,1e-8)%調(diào)用函數(shù)quad求定積分I=3.141592653733437n=61>>[I,n]=quadl(f,0,1,1e-8)%調(diào)用函數(shù)quadl求定積分I=3.141592653589806n=48>>(atan(1)-atan(0))*4%理論值ans=

3.141592653589793>>formatshort第7章MATLAB數(shù)值微分與積分2．自適應(yīng)積分法MATLAB提供了基于全局自適應(yīng)積分算法的integral函數(shù)來求定積分，函數(shù)的調(diào)用格式為：I=integral(filename,a,b)其中，I是計算得到的積分；filename是被積函數(shù)；a和b分別是定積分的下限和上限，積分限可以為無窮大。第7章MATLAB數(shù)值微分與積分例7-5求。先建立被積函數(shù)文件fe.m。functionf=fe(x)f=1./(x.*sqrt(1-log(x).^2));再調(diào)用數(shù)值積分函數(shù)integral求定積分，命令如下：>>I=integral(@fe,1,exp(1))I=1.5708第7章MATLAB數(shù)值微分與積分3．高斯-克朗羅德法MATLAB提供了基于自適應(yīng)高斯-克朗羅德法的quadgk函數(shù)來求振蕩函數(shù)的定積分。該函數(shù)的調(diào)用格式為[I,err]=quadgk(filename,a,b)其中，err返回近似誤差范圍，其他參數(shù)的含義和用法與quad函數(shù)相同。積分上下限可以是?Inf或Inf，也可以是復(fù)數(shù)。如果積分上下限是復(fù)數(shù)，則quadgk在復(fù)平面上求積分。第7章MATLAB數(shù)值微分與積分例7-6求。建立被積函數(shù)文件fsx.m。functionf=fsx(x)f=sin(1./x)./x.^2;調(diào)用函數(shù)quadgk求定積分，命令如下：>>I=quadgk(@fsx,2/pi,+Inf)I=1.0000第7章MATLAB數(shù)值微分與積分4．梯形積分法在科學(xué)實驗和工程應(yīng)用中，函數(shù)關(guān)系表達式往往是不知道的，只有實驗測定的一組樣本點和樣本值，這時，人們就無法使用quad等函數(shù)計算其定積分。在MATLAB中，對由表格形式定義的函數(shù)關(guān)系的求定積分問題用梯形積分函數(shù)trapz，其調(diào)用格式為：I=trapz(X,Y)其中，向量X、Y定義函數(shù)關(guān)系Y

=

f(X)。X、Y是兩個等長的向量：X

=

(x1，x2，…，xn)，Y

=

(y1，y2，…，yn)，并且x1<x2<…<xn，積分區(qū)間是[x1，xn]。第7章MATLAB數(shù)值微分與積分例7-7用trapz函數(shù)計算定積分。命令如下：>>formatlong>>x=0:0.001:1;>>y=4./(1+x.^2);%生成函數(shù)向量>>trapz(x,y)ans=3.141592486923127>>formatshort第7章MATLAB數(shù)值微分與積分5．累計梯形積分在MATLAB中，提供了對數(shù)據(jù)積分逐步累計的函數(shù)cumtrapz。該函數(shù)調(diào)用格式如下。Z=cumtrapz(Y)Z=cumtrapz(X,Y)對于向量Y，Z是一個與Y等長的向量；對于矩陣Y，Z是一個與Y相同大小的矩陣，累計計算Y每列的積分。函數(shù)其他參數(shù)的含義和用法與trapz函數(shù)的相同。例如：>>S=cumtrapz([1:5;2:6]')S=001.50002.50004.00006.00007.500010.500012.000016.0000第7章MATLAB數(shù)值微分與積分7.2.3多重定積分的數(shù)值求解重積分的被積函數(shù)是二元函數(shù)或三元函數(shù)，積分范圍是平面上的一個區(qū)域或空間中的一個區(qū)域。MATLAB提供的integral2、quad2d、dblquad函數(shù)用于求的數(shù)值解，integral3、triplequad函數(shù)用于求的數(shù)值解。函數(shù)的調(diào)用格式為：I=integral2(filename,a,b,c,d)I=quad2d(filename,a,b,c,d)I=dblquad(filename,a,b,c,d,tol)I=integral3(filename,a,b,c,d,e,f)I=triplequad(filename,a,b,c,d,e,f,tol)第7章MATLAB數(shù)值微分與積分例7-8計算二重定積分。建立一個函數(shù)文件fxy.m。functionf=fxy(x,y)globalki;ki=ki+1;%ki用于統(tǒng)計被積函數(shù)的調(diào)用次數(shù)f=exp(-x.^2/2).*sin(x.^2+y);調(diào)用函數(shù)求解，命令如下：>>globalki;>>ki=0;>>I=integral2(@fxy,-2,2,-1,1)%調(diào)用integral2函數(shù)求解I=1.5745>>kiki=22>>ki=0;>>I=quad2d(@fxy,-2,2,-1,1)%調(diào)用quad2d函數(shù)求解I=1.5745>>kiki=20>>ki=0;>>I=dblquad(@fxy,-2,2,-1,1)%調(diào)用dblquad函數(shù)求解I=1.5745>>kiki=1050第7章MATLAB數(shù)值微分與積分例7-9計算三重定積分。命令如下：>>fxyz=@(x,y,z)4*x.*z.*exp(-z.*z.*y-x.*x);%定義被積函數(shù)>>integral3(fxyz,0,pi,0,pi,0,1)ans=1.7328>>triplequad(fxyz,0,pi,0,pi,0,1,1e-7)ans=1.7328第7章MATLAB數(shù)值微分與積分7.3離散傅里葉變換MATLAB提供了一套計算快速傅里葉變換的函數(shù)，它們包括求一維、二維和N維離散傅里葉變換函數(shù)fft、fft2和fftn，還包括求上述各維離散傅里葉變換的逆變換函數(shù)ifft、ifft2和ifftn等。第7章MATLAB數(shù)值微分與積分7.3.1離散傅里葉變換算法簡介在某時間片等距地抽取N個抽樣時間tm處的樣本值f(tm)，且記作f(m)，這里m=0，1，2，…，N-1，稱向量F(k)（k=0，1，2，…，N-1）為f(m)的一個離散傅里葉變換，其中：因為MATLAB不允許有零下標(biāo)，所以將上述公式中m的下標(biāo)均移動1，于是便得相應(yīng)公式：由f(m)求F(k)的過程，稱為求f(m)的離散傅里葉變換，或稱為F(k)為f(m)的離散頻譜。反之，由F(k)逆求f(m)的過程，稱為離散傅里葉逆變換，相應(yīng)的變換公式為：第7章MATLAB數(shù)值微分與積分7.3.2離散傅里葉變換的實現(xiàn)MATLAB提供了對向量或直接對矩陣進行離散傅里葉變換的函數(shù)。下面只介紹一維離散傅里葉變換函數(shù)，其調(diào)用格式為：①fft(X)：返回向量X的離散傅里葉變換。設(shè)X的長度（即元素個數(shù)）為N，若N為2的冪次，則為以2為基數(shù)的快速傅里葉變換，否則為運算速度很慢的非2冪次的算法。對于矩陣X，fft(X)應(yīng)用于矩陣的每一列。②fft(X,N)：計算N點離散傅里葉變換。它限定向量的長度為N，若X的長度小于N，則不足部分補上零；若大于N，則刪去超出N的那些元素。對于矩陣X，它同樣應(yīng)用于矩陣的每一列，只是限定了每一列的長度為N。③fft(X,[],dim)或fft(X,N,dim)：這是對于矩陣而言的函數(shù)調(diào)用格式，前者的功能與fft(X)基本相同，而后者則與fft(X,N)基本相同。只是當(dāng)參數(shù)dim=1時，該函數(shù)作用于X的每一列；當(dāng)dim=2時，則作用于X的每一行。第7章MATLAB數(shù)值微分與積分值得一提的是，當(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ù)或操作方向。第7章MATLAB數(shù)值微分與積分例7-10給定數(shù)學(xué)函數(shù)x(t)=12sin(2π×10t+π/4)+5cos(2π×40t)取N=128，試對t從0~1s采樣，用fft函數(shù)作快速傅里葉變換，繪制相應(yīng)的振幅—頻率圖。在0~1s時間范圍內(nèi)采樣128點，從而可以確定采樣周期和采樣頻率。由于離散傅里葉變換時的下標(biāo)應(yīng)是從0到N-1，故在實際應(yīng)用時下標(biāo)應(yīng)該前移1。又考慮到對離散傅里葉變換來說，其振幅|F(k)|是關(guān)于N/2對稱的，故只須使k從0～N/2即可