臺(tái)灣國(guó)立大學(xué)郭彥甫Matlab教程筆記（17）numerical integration
數(shù)值積分
calculating the numerical value of a definite integral有限整數(shù)

quadrature method （求積分方法）

quadrature method （求積分方法）----approximating the integral by using a finite set of points 通過(guò)有限的點(diǎn)集，近似積分。.quadrature method 把x等分

basic quadrature rules基本求積分方法

1.midpoint rule(zeroth-order approximation)
用矩形去做近似

2.trapezoid rule(first-order approximation)
用梯形去做近似

我們先來(lái)看 Midpoint Rule

對(duì)f(x)做積分，積分區(qū)間是x0~x3，把這個(gè)區(qū)間等分，每個(gè)等分距離是h
取每個(gè)等分小區(qū)間的中點(diǎn)處的函數(shù)值，來(lái)當(dāng)作矩形的高

這個(gè)近似的過(guò)程如下圖

小矩形的面積=底h乘以高f

實(shí)際上怎么用matlab做midpoint rule 呢？

midpoint rule using sum()
舉例求解以下積分：

例程代碼

h=0.05; x=0:h:2; midpoint =(x(1:end-1)+x(2:end))./2;%找所有的中點(diǎn) y=4*midpoint.^3;%被積函數(shù)， s=sum(h*y) %sum求得積分值

來(lái)理解中間這行代碼的目的：
midpoint =(x(1:end-1)+x(2:end))./2;
由于midpoint rule 我們需要把區(qū)間中點(diǎn)值給f（x），所以需要計(jì)算出來(lái)中點(diǎn)值，midpoint,這句代碼就是用來(lái)求中點(diǎn)值的。

執(zhí)行結(jié)果：得到積分值是15.9950

how accurate is it?
積分值是16，用midpoint rule 計(jì)算得到是 15.9950，還算比較精準(zhǔn)

how to improve the accuracy?
步長(zhǎng)h減小

我們?cè)賮?lái)看 trapezoid rule

trapezoid 梯形，不規(guī)則四邊形

和上面的 midpoint rule 近似， 同樣需要等分區(qū)間。這里每個(gè)小區(qū)間近似的函數(shù)值變掉了。
其中，h(f0+f1)/2表示的是梯形面積計(jì)算公式，（上底+下底）*高/2

在matlab中如何使用trapezoid Rule 計(jì)算積分呢？
trapezoid rule using trapz()
同樣的例題

例程代碼：

h=0.05;%等分的小區(qū)間的長(zhǎng)度 x=0:h:2;%等分x y=4*x.^3;%被積函數(shù) s=h*trapz(y)%在小區(qū)間上使用trapz計(jì)算梯形的面積

下圖幫助理解：

計(jì)算結(jié)果：

還有另外一種使用形式

alternative:

h=0.05; x=0:h:2; y=4*x.^3; trapezoid=(y(1:end-1)+y(2:end))/2;%求梯形的（上底+下底）/2 s=h*sum(trapezoid)

再多介紹一種second-order rule :1/3 Simpson’s

這個(gè)近似公式為：（很精準(zhǔn)）（這里是用的兩個(gè)小區(qū)間，所以只有三項(xiàng)，后面會(huì)看到更多的）

積分過(guò)程

Simpson’s rule
計(jì)算過(guò)程：
同樣的例子：

程式碼：

h=0.05; x=0:h:2; y=4*x.^3; s=h/3*(y(1)+2*sum(y(3:2:end-2))+4*sum(y(2:2:end))+y(end))

運(yùn)行結(jié)果：

我們把三種方法做一個(gè)對(duì)比
comparison

用圓圈處的值做逼近

具體區(qū)別如下：

下面講 函數(shù)句柄review of function handles(@)

a handle is a pointer to a function 函數(shù)句柄就是一個(gè)函數(shù)的指針
can be used to pass functions to other functions

通常，一個(gè)function 不能當(dāng)作另外一個(gè)function 的input，函數(shù)句柄可以起到傳遞的作用
舉例子：這里的input 是某一個(gè)function，比如sin

function [y]=xy_plot(input,x) %xy_plot receives the handle of a function and plots that function of x y=input(x); plot(x,y,'r--'); xlabel('x'); ylabel('function(x)'); end

把上面這段代碼儲(chǔ)存成為.m文件，
然后嘗試以下代碼，去呼叫這個(gè)function
xy_plot(sin,0:0.01:2pi);
會(huì)報(bào)錯(cuò)，不能直接給這個(gè)function，需要給這個(gè)function 的指針pointer
正確的用法：
xy_plot(@sin,0:0.01:2pi);
運(yùn)行結(jié)果：我們會(huì)發(fā)現(xiàn)，調(diào)用了這個(gè)sin函數(shù)，繪制了下圖

我們看下一部分numerical integration :integral()

numerical integration on a function from using global adaptive quadrature and default error tolerances

例程：

y=@(x) 1./(x.^3-2*x-5);%讓@(x) 指向被積函數(shù) integral(y,0,2);

代碼的解釋：
integral(y,0,2);%y是被積函數(shù)，需要用函數(shù)句柄來(lái)表示，0和2是積分上下限

這個(gè)積分的結(jié)果：

思考，如果 要算sin（x）在[0,2]的積分怎么做？
同樣的

y=@(x) sin(x);%讓@句柄指向被積函數(shù)sin(x) integral(y,0,2)

double and triple integrals兩重或三重積分

二重積分使用函數(shù)：Integral2()
具體用法：integral(f,pi,2*pi,0,pi)%被積函數(shù)，和二重積分的兩組積分限
計(jì)算下面的兩重積分

例程：

f = @(x,y) y.*sin(x)+x.*cos(y); integral(f,pi,2*pi,0,pi)%被積函數(shù)，和二重積分的兩組積分限

這個(gè)二重積分的計(jì)算結(jié)果：

三重積分用函數(shù)integral3()
具體用法：
integral3(f,0,pi,0,1,-1,1)%參數(shù)是被積函數(shù)，加上三對(duì)積分上下限

計(jì)算下列三重積分：

例程：

f=@(x,y,z) y.*sin(x)+z.*cos(y); integral3(f,0,pi,0,1,-1,1)%參數(shù)是被積函數(shù)，加上三對(duì)積分上下限

上述積分的值為：

【總結(jié)一下】
本文記錄了數(shù)值積分的一些用法。
介紹了三種求數(shù)值積分的方法：用矩形近似sum()函數(shù)，用梯形近似trapz()函數(shù)，還有Simpson方法。然后介紹了函數(shù)作為input的用法， 在此基礎(chǔ)上講解了數(shù)值積分integral()函數(shù)的用法：可以求解一重積分integral()、兩重積分integral2()、三重積分integral3()函數(shù)等。

總結(jié)

以上是生活随笔為你收集整理的台湾国立大学郭彦甫Matlab教程笔记（17）numerical integration的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。