數(shù)值微積分課件

MATLAB數(shù)學建模與仿真數(shù)值微積分2定積分計算的基本公式是牛頓－萊布尼茲公式。但當被積函數(shù)的原函數(shù)不知道時，如何計算？這時就需要利用近似計算。特別是在許多實際應用中，被積函數(shù)甚至沒有解析表達式，而是一條實驗記錄曲線，或一組離散的采樣值，此時只能用近似方法計算定積分。本實驗主要研究定積分的三種近似計算算法：矩形法、梯形法和拋物線法。同時介紹Matlab

計算定積分的相關函數(shù)。

問題背景和實驗目的定積分的近似計算1.極限和連續(xù)數(shù)列極限:>0,N>0，使當n>N時有xn-a<，則函數(shù)極限:如果當xx0時有f(x)A，則連續(xù):如果當xx0時，有f(x)f(x0)

則稱f(x)在x0連續(xù)。閉區(qū)間上連續(xù)函數(shù)必有最大值和最小值。預備知識：微積分當n=0得，微分中值定理

f(x)-f(x0)=f’()(x-x0)

其中是x0與x之間某個值Taylor公式:當f(x)在含有x0某個開區(qū)間內具有直到n+1階的導數(shù)，3.多元函數(shù)微分學

設f(x,y)在點(x0,y0)附近有定義，當(x,y)以任何方式趨向于(x0,y0)時，f(x,y)趨向于一個確定的常數(shù)A，則若A=f(x0,y0),稱f(x,y)在(x0,y0)點連續(xù)f(x,y)在點(x0,y0)的偏導數(shù)分別定義為4.積分

函數(shù)f(x)在區(qū)間[a,b]上的積分定義為其中a=x0<x1<…<xn=b,

xi=xi-xi-1,i(xi-1,xi),i=1,2,…,n若在[a,b]上,F’(x)=f(x),則二重積分定義為8矩形法梯形法拋物線法

數(shù)值積分的常見算法主要內容

Matlab

求積分函數(shù)數(shù)值積分函數(shù)：trapz、quad、dblquad符號積分函數(shù)：int10矩形法n

充分大，x

充分小

通常我們取左點法右點法中點法點可以任意選取，常見的取法有：

左端點,右端點和中點。定積分的近似：12矩形法舉例例：用不同的矩形法計算下面的定積分(取n=100)，

并比較這三種方法的相對誤差。左點法：右點法：中點法：解：h=1/n=0.01,xi=i*h,a=0,b=1,n=100(i=0,1,2,...,100)14定積分幾何意義15

曲邊小梯形的面積可以由直邊小梯形的面積來近似整個曲邊梯形的面積：梯形法16

如果我們n

等分區(qū)間[a,b]，即令：則==>梯形公式梯形法梯形公式與中點公式有什么區(qū)別

?

fuluB.m17解：==>例：用梯形法計算下面定積分(取n=100)，并計算相對誤差梯形法舉例a=0,b=1,n=100,f(x)=1/(1+x2)==>h=1/100=0.01,xi=i*h,yi=f(xi)

相對誤差：182n

等分區(qū)間[a,b]，得用拋物線代替該直線，計算精度是否會更好？

計算每個節(jié)點上的函數(shù)值：拋物線法

在區(qū)間[x0,x2]上，用過以下三點的拋物線來近似原函數(shù)f(x)。20同理可得：相加即得：拋物線法21整理后可得：或辛卜生(Simpson)公式拋物線法公式拋物線法

fuluC.m23矩形法梯形法拋物線法

數(shù)值積分的常見算法Matlab函數(shù)

Matlab

求積分函數(shù)數(shù)值積分函數(shù)：trapz、quad、dblquad符號積分函數(shù)：int24矩形法總結

Matlab

數(shù)值積分函數(shù)：trapz、quad、dblquad梯形法拋物線法26前面的做法例：用梯形法計算下面定積分(取n=100)解：a=0,b=1,n=100,yi

=f(xi)=1/(1+xi2)>>

x=0:1/100:1;>>

y=1./(1+x.^2);>>

trapz(x,y)trapz函數(shù)trapz(x,1./(1+x.^2))trapz舉例27quad(f,a,b,tol)f=f(x)為被積函數(shù)，[a,b]為積分區(qū)間，tol

為計算精度將自變量看成是向量不用自己分割積分區(qū)間可以指定計算精度，若不指定，缺省精度是10-6精度越高，函數(shù)運行的時間越長此處的函數(shù)

f是數(shù)值形式，應該使用數(shù)組運算，即：

.*

./

.\

.^

quad

quad拋物線法28解：>>

quad('1./(1+x.^2)',0,1)>>

quad('1./(1+x.^2)',0,1,1e-10)>>

quad('1./(1+x.^2)',0,1,1e-16)函數(shù)表達式一定要用單引號括起來！涉及的運算一定要用數(shù)組運算！例：用quad

計算定積分：quad舉例30>>

f=inline('4*x*y+3*y^2');>>

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

f

中關于第一自變量的運算是數(shù)組運算，即把x

看成是向量，y

看成是標量。也可以全部采用數(shù)組運算例：計算二重積分>>

dblquad(inline('4*x*y+3*x^2'),-1,1,0,2)>>

dblquad(inline('4*x*y+3*x.^2'),-1,1,0,2)X例：計算二重積分dblquad舉例31例：計算二重積分>>

dblquad(@(x,y)4*x*y+3*x.^2,-1,1,0,2)指定x、y

分別是第一和第二積分變量>>

dblquad(inline('4*x*y+3*x.^2'),-1,1,0,2)被積函數(shù)f(x,y)

的另一種定義方法：匿名函數(shù)>>

dblquad(@(y,x)4*x*y+3*x.^2,-1,1,0,2)下面的命令運行結果和上面的一樣嗎？dblquad舉例32int(f,a,b)

計算

f

關于默認自變量

的定積分，積分區(qū)間為[a,b]。int(f)

計算

f

關于默認自變量

的不定積分。int(f,v,a,b)

計算函數(shù)f

關于自變量v

的定積分，積分區(qū)間為[a,b]int(f,v)

計算函數(shù)

f

關于自變量

v

的不定積分findsym(f,1)int符號積分：int33例：用int

函數(shù)計算定積分：解：>>

symsx;>>

f=1/(1+x^2);>>

int(f,x,0,1)>>

f=sym('1/(1+x^2)');>>

int(f,'x',0,1)>>

int('1/(1+x^2)','x',0,1)或>>

int('1/(1+x^2)',0,1)或或int舉例34double(a)將

a

轉化為雙精度型，若

a

是字符，則取對應的

ASCII碼>>

a=3;>>

double(a)>>

double('a')例：ans=3ans=97相關函數(shù)35梯形法：拋物線法：符號積分法：用Matlab函數(shù)近似計算定積分作業(yè)36拋物線法：符號積分法：作業(yè)用Matlab函數(shù)近似計算二重積分Matlab中幾個數(shù)值積分函數(shù)的比較和優(yōu)缺點

一、Z=trapz(X,Y,dim)

梯形數(shù)值積分，通過已知參數(shù)x,y按dim維使用梯形公式進行積分。例

計算int(sin(x),0,pi)>>x=0:pi/100:2*pi;

>>y=sin(x);

>>z=trapz(x,y)%或者說使用z=pi/100*trapz(y)

z=1.0300e-017>>z=pi/100*trapz(y)二、[q,fcnt]=quad(fun,a,b,tol,trace,p1,p2...)

自適應simpson公式數(shù)值積分，適用于精度要求低，被積函數(shù)平滑性較差的數(shù)值積分。注意事項：

1.被積函數(shù)fun必須是函數(shù)句柄；

2.積分限[a,b]必須是有限的，因此不能為inf；

3.p1為其他需要傳遞的參數(shù)，一般是數(shù)值?？赡芫妫?/p>

1.'Minimumstepsizereached'

意味著子區(qū)間的長度與計算機舍入誤差相當，無法繼續(xù)計算了。原因可能是有不可積的奇點；

2.'Maximumfunctioncountexceeded'

意味著積分遞歸計算超過了10000次。原因可能是有不可積的奇點；

3.'InfiniteorNot-a-Numberfunctionvalueencountered'

意味著在積分計算時，區(qū)間內出現(xiàn)了浮點數(shù)溢出或者被零除。例2計算積分1/(x^3-2*x-p)，其中參數(shù)p=5，積分區(qū)間為[0,2]。>>F=@(x,n)1./(x.^3-2*x-n);

>>Q=quad(@(x)F(x,5),0,2)%或者使用quad(F,0,2,[],[],5)效果是一樣的，%只是前者使用的函數(shù)嵌套。Q=

-0.4605>>quad(F,0,2,[],[],5)ans=

-0.4605三、[q,fcnt]=quadl(fun,a,b,tol,trace,p1,p2...)

自適應Lobatto數(shù)值積分，適用于精度要求高，被積函數(shù)曲線比較平滑的數(shù)值積分。注意事項：

同quad可能警告：

同quad例

計算積分1/(x^3-2*x-p)，其中參數(shù)p=5，積分區(qū)間為[0,2]。

>>F=@(x,p)1./(x.^3-2*x-p);

>>Q=quadl(F,0,2,[],[],5)%或者Q=quadl(@(x)F(x,5),0,2)。Q=

-0.4605四、[q,errbnd]=quadgk(fun,a,b,param1,val1,param2,val2,...)

自適應Gauss-Kronrod數(shù)值積分，適用于高精度和震蕩數(shù)值積分，支持無窮區(qū)間，并且能夠處理端點包含奇點的情況，同時還支持沿著不連續(xù)函數(shù)積分，復數(shù)域線性路徑的圍道積分法。注意事項：

1.積分限[a,b]可以是[-inf,inf]，但必須快速衰減；

2.被積函數(shù)在端點可以有奇點，如果區(qū)間內部有奇點，將以奇點區(qū)間劃分成多個，也就是說奇點只能出現(xiàn)在端點上；

3.被積函數(shù)可以劇烈震蕩；

4.可以計算不連續(xù)積分，此時需要用到'Waypoints'參數(shù)，'Waypoints'中的點必須嚴格單調；

5.可以計算圍道積分，此時需要用到'Waypoints'參數(shù)，并且為復數(shù)，各點之間使用直線連接；

6.param,val為函數(shù)的其它控制參數(shù)，比如上面的'waypoints'就是，具體看幫助。例

計算有奇點積分int(exp(x)*log(x),0,1)。>>F=@(x)exp(x).*log(x);%奇點必須在端點上，否則請先進行區(qū)間劃分。

>>Q=quadgk(F,0,1)Q=

-1.3179例

計算半無限震蕩積分int(x^5*exp(-x)*sin(x),0,inf)。

>>F=@(x)x.^5.*exp(-x).*sin(x);

>>fplot(F,[0,100])%繪圖，看看函數(shù)的圖形

>>[q,errbnd]=quadgk(F,0,inf,'RelTol',1e-8,'AbsTol',1e-12)%積分限中可以有inf，但必須快速收斂q=-15.0000

errbnd=9.4386e-009五、q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)

矩形區(qū)域二重數(shù)值積分，一般區(qū)域二重積分參見NIT(數(shù)值積分工具箱)的quad2dggen函數(shù)。例9計算下面二重積分

>>F=@(x,y)y*sin(x)+x*cos(y);

>Q=dblquad(F,pi,2*pi,0,pi)Q=

-9.8696六q=triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method)

長方體區(qū)域三重數(shù)值積分，注意此時沒有一般區(qū)域的三重積分。例10計算下面三重積分。>>F=@(x,y,z)y*sin(x)+z*cos(x);

>>Q=triplequad(F,0,pi,0,1,-1,1)Q=

2.0000總結(1)quad：采用自適應變步長simpson方法，速度和精度都是最差的，建議不要使用；

(2)quad8：使用8階Newton-Cotes算法，精度和速度均優(yōu)于quad，但在目前版本下已被取消；

(3)quadl：采用lobbato算法，精度和速度均較好，建議全部使用該函數(shù)；

(4)quadg：NIT(數(shù)值積分)工具箱函數(shù)，效率最高，但該工具箱需要另外下載；

(5)quadgk：很有用的函數(shù)，功能在Matlab中最強大；

(6)quad2dggen：一般區(qū)域二重積分，效率很好，需要NIT支持；

(7)dblquad：長方形區(qū)域二重積分；

(8)triplequad：長方體區(qū)域三重積分；

1.導數(shù)、單調性與極值

當f’(x0)>0,函數(shù)在x0點附近是上升的，

f’(x0