第7章 MATLAB數值微分與積分

第7章MATLAB數值微分與積分

7.1數值微分7.2數值積分7.3離散傅里葉變換第7章MATLAB數值微分與積分7.1數值微分7.1.1數值差分與差商任意函數f(x)在x點的導數是通過極限定義的：如果去掉上述等式右端的h→0的極限過程，并引進記號：稱△f(x)、▽f(x)及δf(x)分別為函數在x點處以h（h>0）為步長的向前差分、向后差分和中心差分。第7章MATLAB數值微分與積分當步長h充分小時，稱△f(x)/h、▽f(x)/h及δf(x)/h分別為函數在x點處以h（h>0）為步長的向前差商、向后差商和中心差商。函數f在點x的微分接近于函數在該點的任意種差分，而f在點x的導數接近于函數在該點的任意種差商。第7章MATLAB數值微分與積分7.1.2數值微分的實現有兩種方式計算任意函數f(x)在給定點x的數值導數：第一種方式是用多項式或樣條函數g(x)對f(x)進行逼近（插值或擬合），然后用逼近函數g(x)在點x處的導數作為f(x)在點x處的導數；第二種方式是用f(x)在點x處的某種差商作為其導數。第7章MATLAB數值微分與積分在MATLAB中，沒有直接提供求數值導數的函數，只有計算向前差分的函數diff，其調用格式為：①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時（默認狀態(tài)），按列計算差分；dim=2，按行計算差分。第7章MATLAB數值微分與積分例7-1設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數值微分與積分例7-2設用不同的方法求函數f(x)的數值導數，并在同一個坐標系中做出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求導數dpdpx=polyval(dp,x); %求dp在假設點的函數值dx=diff(f([x,3.01]))/0.01; %直接對f(x)求數值導數gx=g(x); %求函數f的導函數g在假設點的導數plot(x,dpx,x,dx,'.',x,gx,'-') %作圖第7章MATLAB數值微分與積分結果表明，用3種方法求得的數值導數比較接近。第7章MATLAB數值微分與積分7.2數值積分7.2.1數值積分基本原理設：

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

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

=

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

=

(x1，x2，…，xn)，Y

=

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