Matlab基礎及其應用教程07課件

1、第7章數(shù)值微積分與常微分方程求解【本章學習目標】掌握數(shù)值微分的實現(xiàn)方法。掌握數(shù)值積分基本原理和實現(xiàn)方法。了解離散傅里葉變換的原理及實現(xiàn)方法。掌握常微分方程的數(shù)值求解方法。7.1 數(shù)值微分稱f(x)、f(x)及f(x)分別為函數(shù)在x點處以h（h0）為步長的向前差分、向后差分和中心差分。稱f(x)/h、f(x)/h及f(x)/h分別為函數(shù)在x點處以h（h0）為步長的向前差商、向后差商和中心差商。7.1.1 數(shù)值差分與差商c1多項式求導法用多項式或樣條函數(shù)g(x)對f(x)進行逼近（插值或擬合），然后用逼近函數(shù)g(x)在點x處的導數(shù)作為f(x)在點x處的導數(shù)。2用diff函數(shù)計算差分用f(x)在點x

2、處的某種差商作為其導數(shù)。在MATLAB中，計算向前差分的函數(shù)diff，調用格式如下。DX = diff(X)：計算向量X的向前差分。DX = diff(X,n)：計算向量X的n階向前差分。DX = diff(X,n,dim)：計算矩陣X的n階差分，dim = 1時（默認狀態(tài)），按列計算差分；dim = 2，按行計算差分。7.1.2 數(shù)值微分的實現(xiàn)7.1 數(shù)值微分【例7.1】設f(x)=sinx，用不同的方法求函數(shù)f(x)的數(shù)值導數(shù)，并在同一個坐標系中繪制f(x)的三種方法所得導數(shù)曲線。x=0:pi/24:pi;%用5次多項式p擬合f(x),并對擬合多項式p求導數(shù)dp在假設點的函數(shù)值p=poly

3、fit(x,sin(x),5);dp=polyder(p);dpx=polyval(dp,x);%直接對sin(x)求數(shù)值導數(shù)dx=diff(sin(x,pi+pi/24)/(pi/24);%求函數(shù)f的導函數(shù)g在假設點的導數(shù)gx=cos(x); plot(x,dpx,b-,x,dx,ko,x,gx,r+);7.1 數(shù)值微分【例7.2】生成一個5階魔方矩陣，按列進行差分運算。M=magic(5)M= 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9DM=diff(M) %計算V的一階差分DM= 6 19 6 6 1

4、 19 1 6 6 6 6 6 6 1 19 1 6 6 19 67.2 數(shù)值積分基本思想是將整個積分區(qū)間a，b分成n個子區(qū)間xi，xi+1，i=1，2，n，其中x1=a，xn+1=b。這樣求定積分問題就分解為下面的求和問題：矩形法是用矩形面積近似曲邊梯形的面積，如圖7.2（a）所示；梯形法是用斜邊梯形面積近似曲邊梯形的面積，如圖7.2（b）所示；而辛普生法是用拋物線近似曲邊。7.2.1 數(shù)值積分的原理7.2 數(shù)值積分1梯形法2自適應辛普生法3高斯-克朗羅德（Gauss-Kronrod）法7.2.1 數(shù)值積分的原理7.2 數(shù)值積分1自適應積分算法基于全局自適應積分算法的integral函數(shù)來求

5、定積分。函數(shù)的調用格式為q = integral(fun,xmin,xmax) q = integral(fun,xmin,xmax,Name,Value)fun是被積函數(shù)，xmin和xmax分別是定積分的下限和上限。選項用于設置積分器的屬性7.2.2 定積分的數(shù)值求解實現(xiàn)7.2 數(shù)值積分【例7.3】求（1）建立被積函數(shù)文件fe.m。function f=fe(x)f=x.*sin(x)./(1+abs(cos(x);（2）調用數(shù)值積分函數(shù)integral求定積分。 q=integral(fe,0,pi)q = 2.1776例7.3也可使用以下命令求解：q=integral(x)( x.*si

6、n(x)./(1+abs(cos(x),0,pi)7.2 數(shù)值積分 a=8755; b=6810; format long funLength=(x)sqrt(a2.*sin(x).2+ b2.*cos(x).2); L=4*integral(funLength, 0, pi/2)L = 4.908996526868900e+047.2 數(shù)值積分2高斯-克朗羅德法quadgk函數(shù)求振蕩函數(shù)的定積分。該函數(shù)的調用格式為q = quadgk(fun, xmin,xmax)q,errbnd = quadgk(fun, xmin,xmax,Name,Value)errbnd返回近似誤差范圍。積分上下限

7、可以是Inf或Inf，如果積分上下限是復數(shù)，則quadgk在復平面上求積分。7.2.2 定積分的數(shù)值求解實現(xiàn)7.2 數(shù)值積分【例7.5】求定積分 format long; I=quadgk(x)exp(x).*log(x),0,1)I= 1.3179021624140817.2 數(shù)值積分3梯形積分法函數(shù)trapz對的離散數(shù)據(jù)用梯形法求定積分，調用格式如下。（1）T= trapz(Y)用于求均勻間距的積分。通常，Y是向量，采用單位間距（即間距為1），計算Y的近似積分。若Y是矩陣，則輸出參數(shù)T是一個行向量，T的每個元素分別存儲Y的每一列的積分結果。（2）T=trapz(X, Y)用于求非均勻間距的

8、積分。X、Y滿足函數(shù)關系Y=f(X)， 按X指定的數(shù)據(jù)點間距，對Y求積分。7.2.2 定積分的數(shù)值求解實現(xiàn)7.2 數(shù)值積分【例7.6】從地面發(fā)射一枚火箭，表7.2記錄了080秒火箭的加速度。試求火箭在第80秒時的速度。 t=0:10:80; a=30.00,31.63,33.44,35.47,37.75,40.33,43.29,46.69,50.67; v=trapz(t,a)v = 3.0893e+037.2 數(shù)值積分4累計梯形積分在MATLAB中，提供了對數(shù)據(jù)積分逐步累計的函數(shù)cumtrapz。該函數(shù)調用格式如下。（1）Z = cumtrapz(Y)用于求均勻間距的累計積分。（2）Z =

9、cumtrapz(X,Y)用于求非均勻間距的累計積分。7.2.2 定積分的數(shù)值求解實現(xiàn)7.2 數(shù)值積分【例7.6】從地面發(fā)射一枚火箭，表7.2記錄了080秒火箭的加速度。試求火箭在第80秒時的速度。計算例7.6中各時間點的速度，命令如下： v=cumtrapz(t,a) v=cumtrapz(t,a)v = 1.0e+03 * 1 至 7 列 0 0.3081 0.6335 0.9780 1.3442 1.7346 2.1526 8 至 9 列 2.6026 3.08947.2 數(shù)值積分integral2、quad2d函數(shù)用于求二重積分的數(shù)值解，integral3函數(shù)用于求三重積分的數(shù)值解。函

10、數(shù)的調用格式為：q = integral2(fun,xmin,xmax,ymin,ymax,Name,Value)q = quad2d(fun, xmin,xmax,ymin,ymax,)q,errbnd = quad2d(fun, xmin,xmax,ymin,ymax, Name,Value)q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax)q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax,Name,Value)fun為被積函數(shù)，xmin, xmax為x的積分區(qū)域，ymin, ymax為y的積分區(qū)域

11、，zmin, zmax為z的積分區(qū)域，選項Name的用法及可取值與函數(shù)integral相同。輸出參數(shù)q返回積分結果，errbnd用于返回計算誤差。7.2.3 多重定積分的數(shù)值求解實現(xiàn)7.2 數(shù)值積分【例7.7】計算二重定積分 fxy = (x,y) exp(-x.2/2).*sin(x.2+y); I= integral2(fxy,-2,2,-1,1)I = 1.5745【例7.8】計算三重定積分 fxyz = (x,y,z)4*x.*z.*exp(-z.*z.*y-x.*x); I3 = integral3(fxyz,0,pi,0,pi,0,1)I3= 1.73287.3 離散傅立葉變換ff

12、t函數(shù)用于求離散傅立葉變換，調用格式如下。fft(X,n,dim)通常，X是向量，若X是矩陣，fft(X)應用于矩陣的每一列。選項n指定變換的數(shù)據(jù)點數(shù)，若X的長度小于n，則不足部分補上零；若大于n，則刪去超出n的那些元素。選項dim指定操作方向，dim=1時，該函數(shù)作用于X的每一列；dim=2時，則作用于X的每一行。7.2.2 定積分的數(shù)值求解實現(xiàn)7.3 離散傅立葉變換【例7.9】給定數(shù)學函數(shù)x(t)= 3cos(220t+/3)+5sin(250t)，在01s時間范圍內采樣128點，用fft作快速傅立葉變換，繪制相應的振幅頻率圖。N=128; % 采樣點數(shù)T=1; % 采樣時間終點t=lin

13、space(0,T,N); x=3*cos(2*pi*20*t+pi/3)+5*sin(2*pi*50*t); % 求各采樣點樣本值dt=t(2)-t(1); % 采樣周期fs=1/dt; % 采樣頻率(Hz)X=fft(x); % 計算x的快速傅立葉變換F=X(1:N); f=fs/N*(0:N-1); % 得到信號的頻率成分plot(f,abs(F),-*) % 繪制振幅頻率圖xlabel(Frequency);ylabel(|F(k)|)7.4 常微分方程的數(shù)值求解4階龍格-庫塔（Runge-Kutta）公式為7.4.1 龍格-庫塔法簡介7.4 常微分方程的數(shù)值求解1.求解函數(shù)求解常微分

14、方程數(shù)值解的函數(shù)用法相同，其基本調用格式為T,Y = solver(odefun,tspan,y0,options)T,Y,TE,YE,IE = solver(odefun,tspan,y0,options)其中，solver是根據(jù)待求解問題的性質選用的求解函數(shù)。7.4.2 常微分方程數(shù)值求解的實現(xiàn)7.4 常微分方程的數(shù)值求解2.求解函數(shù)的參數(shù)（1）輸出參數(shù)T為列向量，返回時間點；Y為矩陣，Y的每一行向量返回T中各時間點相應的狀態(tài)。TE為事件發(fā)生的時刻，YE為對應時刻的解，IE為觸發(fā)事件的索引。（2）輸入?yún)?shù)odefun為函數(shù)文件名或匿名函數(shù)，函數(shù)文件頭通常采用f(t，y)的形式，即形參t為時

15、間參量，形參y為待求解問題的自變量。tspan指定求解區(qū)間，用二元向量t0 tf表示。y0是初始狀態(tài)列向量。options用于設置積分求解過程和結果的屬性。7.4.2 常微分方程數(shù)值求解的實現(xiàn)7.4 常微分方程的數(shù)值求解【例7.10】設有初值問題試求其數(shù)值解，并與精確解相比較，精確解為（1）建立函數(shù)文件funst.m。function yp=funst(x,y)yp= y.*tan(x) + sec(x);（2）求解微分方程。 y0=pi/2; x,y=ode23(funst,0,1,y0);%求數(shù)值解 y1=(x+pi/2)./cos(x); %求精確解 x,y,y1ans= 0 1.570

16、8 1.5708 0.1000 1.6792 1.6792 0.2000 1.8068 1.8068 0.3000 1.9583 1.9583 0.4000 2.1397 2.1397 0.5000 2.3596 2.3597 0.6000 2.6301 2.6302 0.7000 2.9689 2.9690 0.8000 3.4027 3.4029 0.9000 3.9745 3.9748 1.0000 4.7573 4.75817.4 常微分方程的數(shù)值求解【例7.11】求描述振蕩器的Van der Pol方程：令x1=y，x2=y，則可寫出Van der Pol方程的狀態(tài)方程：（1）建立函數(shù)文件verderpol.m。function xprime=verderpol(t,x)global mu;xprime=x(2); mu*(1-x(1)2)*x(2)-x(1);（2）求解微分方程。 global mu; mu=2; y0=1;0; t,x=ode45(verderpol,0,20,y0); subplot(1,2,1);plot(t,x); %繪制系統(tǒng)時間響應曲線 subplot(1,2,2);plot(x(:,1),x(:,2) %繪制系統(tǒng)相平面曲線7.4 常微分方程的數(shù)值求