Matlab課程設計--數(shù)值微積分

1、；.數(shù)值微積分引言 在微分中，函數(shù)的導數(shù)是用極限來定義的，如果一個函數(shù)是以數(shù)值給出的離散形式，那么它的導數(shù)就無法用極限運算方法求得，當然也就更無法用求道方法去計算函數(shù)在某點處的導數(shù)。 一般來說，函數(shù)的導數(shù)依然是一個函數(shù)。設函數(shù)f(x)的導數(shù)f(x)=g(x),高等數(shù)學關心的是g（x）的形式和性質，而數(shù)值分析關心的問題是怎樣的計算g（x）在一串離散點X=（x1,x2,xn）的近似值G=（g1,g2,.gn）以及所計算的近似值有多大的誤差。1.    數(shù)值差分與差商任意函數(shù)f(x)在x點的導數(shù)是通過極限定義的： f(x)=lim f(x+h)-f(x)/h f(x)=

2、lim f(x)-f(x-h)/h f(x)=lim f(x+h/2)-f(x-h/2)/h 上述式子中，均假設h0，如果去掉上述等式右端的h0的極限過程，并引進記號：f(x)=f(x+h)-f(x)f(x)=f(x)-f(x-h)f(x)=f(x+h/2)-f(x-h/2) 稱f(x) f(x)、f(x)分別為函數(shù)在x點的以h（h>0）為步長的向前差分、向后差分、向中差分。當步長h充分的小時，有  f(x）f(x)/h f(x）f(x)/h f(x）f(x)/h 和差分一樣，稱f(x)/h, f(x)/h,、f(x)/h分別為函數(shù)在x點一h（h>0）

3、為步長的向前差商、向后差商、向中差商。當步長h充分的小時，函數(shù)f在x點微分接近于函數(shù)在該點的任意種差分，而f在點x的導數(shù)接近于函數(shù)在該點的任意種差商。 2.    數(shù)值微分的實現(xiàn)有兩種方式計算任意函數(shù)f（x）在給定點x點數(shù)值導數(shù)。第一種方式是用多項式或樣條函數(shù)g（x）對f(x)進行逼近，然后用逼近函數(shù)g(x)在點x的導數(shù)作為f(x)在點x的導數(shù)。第二種方式是用f（x）在點x處的某種差商作為其導數(shù)。在MATLAB中，沒有直接提供數(shù)值導數(shù)的函數(shù)，只有計算向前差分的函數(shù)diff，其調用格式為： Diff(X),輸入?yún)?shù)x可為矢量或矩陣。若X為矢量，則返回X(2

4、)-X(1),X(3)-X(2),X(n)-X(n-1);若X為矩陣，則返回矩陣每行的差分，即X(2:m，：)-X(1:m-1,:).通常，diff(X)返回沿X的第一個成對維的差值。 Y=diff(X,n):遞歸調用diff函數(shù)n次，生成n次差分。 Y=diff(X,n,dim): 返回X中第dim維的n次差分或導數(shù)值 設計1如下： X=3 7 5;0 4 2; Diff(x) ans= -3 -3 -3 設計2如下： x=3 7 5;0 4 2; Diff=(x,2) ans=0       

5、0;         0   設計3如下：x=3 7 5；0 4 2；diff=(x,2,1)ans= empty matrix: 0-by-3  設： 用不同的方法求函數(shù)f(x)的數(shù)值導數(shù)，并在同一個坐標系中畫出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、/6)+5); x=-3:0.01:3; p=polyfit(x,f(x),5); 用5次多項式p擬合f(x) dp=polyder(p); 對擬合多項式p求導數(shù)dpdpx=polyval(dp,x); 求dp在假設點的函數(shù)值dx=diff(f(x,3.01)/0.01; 直接對f(x)求數(shù)值導數(shù) gx=g(x); 求函數(shù)f的導函數(shù)g在假設點的導數(shù)plot(x,dpx,x,dx,.,x,gx,-); 作圖  程序運行后得到圖6.3所示的圖形。結果表明，用3種方法求得的數(shù)值導數(shù)比較近似。        &#

7、160;  MATLAB語言提供了計算給定向量差分的函數(shù)diff(),其調用的方法是dy=diff(y).假設向量y是yi,i=1,2,3,n構成，則經(jīng)過diff()函數(shù)處理后將得出一個新的向量：yi+1-yi,i=1,2,3,n,顯然新得出的向量長度比原來向量的長度y要小1.如：  >>v=vander(1:6) v= 1 1 1 1 1 1 32 16 8 4 2 1 243 81 27 9 3 1 1024 256 64 16 4 1 3125 625 125 25 5 1 7776 1296 216 36 6 1  >>d

8、iff(v) ans= 31 15 7 3 1 0 211 65 19 5 1 0 781 175 37 7 1 0 2101 369 61 9 1 0 4651 671 91 11 1 0 可見，diff（）函數(shù)對矩陣的每一列都進行了差分運算，故而結果矩陣的列數(shù)是不變的，只是有行數(shù)減1.  同時，在MATLAB中，還提供了gradient()函數(shù)的調運，這函數(shù)可以直接用來求一個矩陣的二維差分，該函數(shù)的調用格式為：  dx,dy=gradient(A) >> v=vander(1,6),dx,dy=gradieng(v) v= 1 1 1 1 1 1

9、32 16 8 4 2 1 243 81 27 9 3 1 1024 256 64 16 4 1 3125 625 125 25 5 1 7776 1296 216 36 6 1 dx= 1.0e=003* 0 0 0 0 0 0 -0.0160 -0.0120 -0.0060 -0.0030 -0.0015 -0.0010 -0.1620 -0.1080 -0.0360 -0.0120 -0.0040 -0.0020 -0.7600 -0.4800 -0.1200 -0.0300 -0.0075 -0.0030 -2.5000 -1.5000 -0.3000 -0.0600 -0.0120

10、-0.0040 -6.4800 -3.7800 -0.6300 -0.1050 -0.0175 -0.0050 dy= 31 15 7 3 1 0 121 40 13 4 1 0 496 120 28 6 1 0 1441 272 49 8 1 0 3376 520 76 10 1 0 4651 671 91 11 1 0   6.2.2 數(shù)值積分 在科學實驗和生產中，經(jīng)常要求函數(shù)f(x)在區(qū)間a,b上的定積分： 在高等數(shù)學中，計算積分依靠微積分基本定理，只要找到被積函數(shù)f(x)的原函數(shù)F（x），則可以用牛頓萊布尼茨公式： 來求出定積分。但是，在有些情況

11、下，應用牛頓萊布尼茨公式往往有困難，例如，當被積函數(shù)的原函數(shù)無法用初等函數(shù)表示或被積函數(shù)為僅知離散點處函數(shù)值的離散函數(shù)時。 1.    數(shù)值積分的基本原理數(shù)值積分研究定積分的數(shù)值求解方法。設 I1= 高等數(shù)學中知道，當f(x)-p(x) <時，I1-I2<(b-a).這說明。當充分的小時，可用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上的積分值。一

12、般選擇被積分函數(shù)的選擇一次多項式時，稱為梯形公式。選擇二次多項式時，為辛譜森（simpson）公式。在MATLAB中，我們通常可以使用梯形求積，Smpson求積，Lobotta求積以及二重積分和三重積分運算，另外還有Gass求積和Romberg求積算法。這些方法的基本思想就是把整個的積分空間分割成為若干個子空間，每個空間上的函數(shù)積分可以求取，因而在整個的空間函數(shù)上積分可以得到，MATLAB基于這種思想采用自適應變步長的方法給出quad()函數(shù)來求積分?，F(xiàn)在分別舉例說明：1.    梯形求積 用trapz函數(shù)進行梯形數(shù)值積分。 設：   的精確值為2，下面

13、用trapz函數(shù)在均勻間隔的網(wǎng)絡上求取該積分的數(shù)值近似。 X=0:pi/100:pi; Y=sin(X); Z=trapz(X,Y) Z= 1.9998 對于非均勻的間隔情況，例如： X=sort(rand(1.101)*pi); Y=sin(X); Z=trapz(X,Y) Z=1.9981  2.Simpson求積用quad和quad8函數(shù)進行自適應遞歸Simpson求積。要求積分具有下面的形式： Q=以quad函數(shù)為例，quad和quad8函數(shù)具有下面的語法形式。         

14、;     q=quad(fun,a,b);求函數(shù)fun從a到b的積分近似，誤差為le-6.fun為M文件函數(shù)或匿名函數(shù)的句柄，函數(shù)x=fun(x)應該接受矢量變量x并返回結果矢量y。              q=quad(fun,a,b):使用絕對誤差容限tol代替默認的容限1.0e-6。Tol的值更大時，需要的計算次數(shù)少，計算快。       

15、0;      q=quad(fun,a,b,tol,trace):在遞歸過程中使用非0的trace參數(shù)。              q.fcnt=quad(): 返回函數(shù)計算次數(shù) 用quad函數(shù)計算下面積分的數(shù)值積分 首先編寫函數(shù)的M文件myfun.m。 Function y=myfun(x) Y=1./(x.3-2*x-5); 然后在命令窗口鍵入 Q=quadruple（myfun.0,2） Q= -0.4605

16、 也可以直接使用匿名函數(shù)形式，如： F=（x）1.9(x.3-2*x-5)； Q=quad(F,0,2); 用以上兩種方法求：  先建立一個函數(shù)文件ex.m: Function ex=ex(x) ex=exp(-x.2);然后在MATLAB命令窗口，輸入命令： Format long I=quad(ex,0,1) I= 0.74682418072642 I=quadl(ex,0,1) I= 0.74682413398845 也可以不建立關于被積函數(shù)文件，而使用語句函數(shù)求解，命令如下： q=inline(exp(-x.2); I=quadl(g,0,1) I= 0.7468241339

17、8845 Format short 讀者可以用解析方法計算下本例，并可以比較一下。同樣。對于下列函數(shù) Humps(x)=   可編寫名為humps.m的M文件，內容是 Functions y=humps(x) y=1./(x-.3).2+0.01)+1./(x-.9).2+0.04)-6； 該函數(shù)的圖形如下： 生成的方法如下： X:-1:0.01:2;Polt(x,humps(x)                現(xiàn)在

18、需要求humps從0到1積分，可使用下面的命令： >>q=quad(humps,0,1) q= 29.8583 注意quad的一般調用格式為 y,n=quad(F,a,b,tol)其中F為描述被積分的字符變量，一般為一個F.m函數(shù)文件名，該函數(shù)的一般格式為y=F(x);a,b為上下限；tol為變步長積分用的誤差限。 MATLAB中常見的一元函數(shù)數(shù)值積分指令如下 表2  quad自適應辛譜森積分Quad8牛頓8段積分quadl自適應Lobtto積分trapz梯形法求積分sum等寬矩形法求積分fnint樣條函數(shù)求積分 其中quad()與quad8()這兩者的調用格式

19、完全一致，只不過在函數(shù)quad8（）中tol的默認值較低。該算法可以更加精確的求出積分的值，但現(xiàn)在在MATLAB一般用新引進的函數(shù)quadl（）來代替quad8(). 現(xiàn)在在來介紹下其他的方法3. LOBATTO求積 用quadl函數(shù)進行Lobatto求積，該函數(shù)的語法格式為：              q=quad(fun,a,b);求函數(shù)fun從a到b的積分近似，誤差為le-6.fun為M文件函數(shù)或匿名函數(shù)的句柄，函數(shù)x=fun(x)應該接受矢量變量x并返回結果

20、矢量y。              q=quad(fun,a,b):使用絕對誤差容限tol代替默認的容限1.0e-6。Tol的值更大時，需要的計算次數(shù)少，計算快。              q=quad(fun,a,b,tol,trace):在遞歸過程中使用非0的trace參數(shù)。     

21、         q.fcnt=quad(): 返回函數(shù)計算次數(shù) 例： 用quadl函數(shù)求下面積分的數(shù)值積分。    首先編寫函數(shù)的M文件myfun.m。 Function y=myfun(x) Y=1./(x.3-2*x-5); 然后在命令窗口鍵入 Q=quadruple（myfun.0,2） Q= -0.4605 也可以直接使用匿名函數(shù)形式，如： F=（x）1.9(x.3-2*x-5)； Q=quad(F,0,2)  二重積分用dblquad函數(shù)對二重積分進行數(shù)值計算，該

22、函數(shù)的語法格式為：              q=dblquad(fun,xmin,xmax,ymin,ymax);調用quad函數(shù)對fun(x,y)函數(shù)進行二重積分運算積分區(qū)間為xmianxxmax;yminyymax。Fun參數(shù)是一個M文件函數(shù)或匿名函數(shù)的句柄。Fun(x,y)函數(shù)必須接受矢量x與標量y，并返回積分值組成的矢量。           &#

23、160;  q=dblquad(fun,xmin,xmax,ymin,ymax,tol):用容限tol代替默認的1.0e-6.              Q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method):用method參數(shù)指定積分方法，代替默認的quad函數(shù)。Method還可以指定為quadl或用戶指定的函數(shù)句柄，該函數(shù)必須與quad和quadl具有相同的調用格式。考慮下面的二重積分問題：  使用MATL

24、AB提供的dblquad函數(shù)就直接可以求出上述積分的數(shù)值解，該函數(shù)的調用格式為：I=dblquad(f,a,b,c,d,tol.trace)該函數(shù)求f(x,y)在a,b×c,d區(qū)域上的二重積分。 例：用dblquad函數(shù)求下面二重積分的數(shù)值近似。   首先建立函數(shù)的M文件integrnd.m,代碼為： Function z=integrnd(x,y) z=y*sin(x)+x*cos(y); 然后在命令窗口鍵入 Q=dblquad(integend,pi,2*pi,0,pi) Q= -9.8696   例： 計算二次定積分 （1）  

25、  建立一個函數(shù)文件fxy.mFunction f=fxy(x,y)Global ki;Ki=ki+I;F=exp(-x.2/2).*sin(x.2+y);(2)      調用dblquad函數(shù)求解Global ki;ki=0;I=dblquad(fxy,-2,2,-1,1)KiI= 1.57449318974494Ki= 1038如果使用inline函數(shù)，則命令如下： f=inline(exp(-x.2/2).*sin(x.2+y),x,y); I=dblquad(f,-2,2m-1,1) I= 1.57449318974494  &