臺灣國立大學郭彥甫Matlab教程筆記（17）numerical integration
數值積分
calculating the numerical value of a definite integral有限整數

quadrature method （求積分方法）

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

basic quadrature rules基本求積分方法

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

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

我們先來看 Midpoint Rule

對f(x)做積分，積分區間是x0~x3，把這個區間等分，每個等分距離是h
取每個等分小區間的中點處的函數值，來當作矩形的高

這個近似的過程如下圖

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

實際上怎么用matlab做midpoint rule 呢？

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

例程代碼

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

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

執行結果：得到積分值是15.9950

how accurate is it?
積分值是16，用midpoint rule 計算得到是 15.9950，還算比較精準

how to improve the accuracy?
步長h減小

我們再來看 trapezoid rule

trapezoid 梯形，不規則四邊形

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

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

例程代碼：

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

下圖幫助理解：

計算結果：

還有另外一種使用形式

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

這個近似公式為：（很精準）（這里是用的兩個小區間，所以只有三項，后面會看到更多的）

積分過程

Simpson’s rule
計算過程：
同樣的例子：

程式碼：

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))

運行結果：

我們把三種方法做一個對比
comparison

用圓圈處的值做逼近

具體區別如下：

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

a handle is a pointer to a function 函數句柄就是一個函數的指針
can be used to pass functions to other functions

通常，一個function 不能當作另外一個function 的input，函數句柄可以起到傳遞的作用
舉例子：這里的input 是某一個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

把上面這段代碼儲存成為.m文件，
然后嘗試以下代碼，去呼叫這個function
xy_plot(sin,0:0.01:2pi);
會報錯，不能直接給這個function，需要給這個function 的指針pointer
正確的用法：
xy_plot(@sin,0:0.01:2pi);
運行結果：我們會發現，調用了這個sin函數，繪制了下圖

我們看下一部分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) 指向被積函數 integral(y,0,2);

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

這個積分的結果：

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

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

double and triple integrals兩重或三重積分

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

例程：

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

這個二重積分的計算結果：

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

計算下列三重積分：

例程：

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

上述積分的值為：

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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。