何俊俊，蘇岐芳

（臺州學院 數(shù)學與信息工程學院，浙江 臨海 317000）

數(shù)值積分的迭代方法及應(yīng)用

何俊俊，蘇岐芳＊

（臺州學院 數(shù)學與信息工程學院，浙江 臨海 317000）

在工程和科學計算中，經(jīng)常會遇到各種類型的積分問題.對于被積函數(shù)過于復雜，其原函數(shù)很難求得，甚至原函數(shù)根本就不是初等函數(shù)；或不知道被積函數(shù)的解析式，而只給出被積函數(shù)在有限個點處的函數(shù)值等情況，需要利用數(shù)值積分方法求積分的近似值．給出了兩種逐次分半求積算法和二重積分的復合梯形算法，并利用這些方法解決了幾類實際問題.

數(shù)值積分；算法；收斂速度；MATLAB程序

1 引言

數(shù)值積分在金融數(shù)學，計算機圖形學，積分方程，工程計算等許多科學領(lǐng)域都有非常重要的應(yīng)用，科學和工程計算中的許多實際問題最終也都會歸結(jié)為數(shù)值積分的計算問題［1］.

由微積分定理，若被積函數(shù)f（x）在區(qū)間［a，b］上連續(xù)，只要能找到f（x）的一個原函數(shù)F（x），便可以利用牛頓-萊布尼茨公式求得積分值.然而，在實際應(yīng)用中常常會碰到一些困難，導致不能用牛頓-萊布尼茨公式直接求得定積分的值．其中，可能出現(xiàn)的情況有［1-2］：

（1）某些函數(shù)，其原函數(shù)不能用初等函數(shù)表示；

（2）函數(shù)f（x）的結(jié)構(gòu)復雜，其原函數(shù)很難求得；

（3）函數(shù)f（x）的表達式不明顯，用表或圖的功能，只給出苦干離散點上的函數(shù)值.

一般地，可根據(jù)給定的節(jié)點構(gòu)造機械求積公式：

其中xk是求積節(jié)點；Ak是求積系數(shù)，也稱伴隨節(jié)點xk的權(quán).

2 復合求積公式

通常把積分區(qū)間等分成若干子區(qū)間，在每個小區(qū)間上利用牛頓-柯特斯公式，再把每個小區(qū)間上的結(jié)果相加作為整個區(qū)間上的積分近似值．本文考慮逐次分半迭代方法．

2.1 梯形遞推公式

將積分區(qū)間［a，b］分為n等份，則有復合梯形公式：

將每個小區(qū)間再次二等分后，得梯形遞推公式［2，8］

2.2 辛普森遞推公式

將積分區(qū)間［a，b］分為n等份，則有復合辛普森公式：

將每個小區(qū)間再次二等分后，

即得辛普森遞推公式

2.3 二重積分求積公式

其中aij是下面矩陣A的對應(yīng)元素：

3 求積迭代算法

算法1：逐次分半梯形法

（1）取k=0，h=b-a，計算T1；

（2）計算第一次分半梯形值T2；

（3）計算誤差err=T2-T1；

（4）如果err＜ε（預(yù)先給定的精度），則停止；否則令k+1→k，T2→T1；

（5）計算T2k，并令T2k→T2轉(zhuǎn)步驟（3）．

算法2：逐次分半辛普森法

（1）取k=0，h=b-a計算S1；

（2）計算第一次分半辛普森值S2；

（3）計算誤差err=S2-S1；

（4）如果err＜ε（預(yù)先給定的精度），則停止；否則令k+1→k，S2→S1；

（5）計算S2k，并令S2k→S2轉(zhuǎn)步驟（3）．

4 數(shù)值積分的應(yīng)用

4.1 衛(wèi)星軌道問題

我國第一顆人造衛(wèi)星近地點距離h=439km，遠地點距離H=2384km，地球衛(wèi)星軌道是一個橢圓，試求衛(wèi)星軌道的周長（地球半徑為R=6371km），且軌道周長公式為

由問題可知，軌道的周長與衛(wèi)星軌道的半軸長a以及地球中心與軌道中心（橢圓中心）的距離c有關(guān)，地球半徑為R=6371km，因而a，c可由近地點、遠地點、地球半徑求出：

在MATLAB-R2011a環(huán)境下，分別用逐次分半梯形法、逐次分半辛普森法和龍貝格算法［2，8］求軌道周長的近似值，要求誤差小于5×10-8，并比較各種算法所需分半次數(shù)及所需CPU時間．運行結(jié)果如表1．

表1 三種算法的計算結(jié)果Table.1 Te results of three algorithm

由輸出結(jié)果看到，對于給定的數(shù)據(jù)，三種算法收斂速度都很快，所需CPU時間基本相同．

4.2 男大學生身高問題

據(jù)資料報道,中國大學生男子組的平均身高是170cm，高度在150cm組至190cm組之間的人大致有99.7%．若是把［150cm，190cm］平均分成20個高度區(qū)間，那男子組身高在每一高度區(qū)間的分布情況如何呢？并且，你能求出身高中等（165cm至175cm之間）的人占男子組的百分比會超過60%嗎？

由問題可知，身高受很多因素影響，通常它是一個服從正態(tài)分布N（u，σ）的隨機變量。正態(tài)分布的概率密度函數(shù)φ（x）為：

則密度函數(shù)φ（x）在區(qū)間［a，b］上的定積分的值能夠較好的代表大學生男子組的身高在區(qū)間［acm，bcm］的分布，在密度函數(shù)φ（x）中的兩個參數(shù)μ、σ分別是正態(tài)分布中的均值與標準差.

據(jù)了解，中國大學生男子組的平均身高是170cm，那么正態(tài)分布選取的均值u=170cm，而“該男子組中約有99.7%的人身高在150cm組至190cm組之間”符合了正態(tài)分布N（u，σ）的“3σ規(guī)則”的條件，即P （μ-3σ＜x≤μ+3σ）=99.7%，得到

將［150cm，190cm］等分成20個高度小區(qū)間，即

在各個小區(qū)間內(nèi)的身高分布情況對應(yīng)于積分值

身高在165cm至175cm之間的人占該男子組的百分比為

即身高中等的人占男子組的百分比．

身高分布函數(shù)φ（x）的圖像如圖1．現(xiàn)用數(shù)值積分方法求積分（8）的近似值，要求誤差小于10-4．

（1）逐次分半梯形法計算結(jié)果，見表2.

表2 逐次分半梯形法計算結(jié)果Table.2 Results by using composite trapezoidal formula

其中n為分半次數(shù)．由（9）式計算得

即身高中等（165cm，175cm）的大學生約為54.67%，不足60%．梯形求積公式計算結(jié)果的分布曲線見圖2．

（2）逐次分半辛普森算法計算結(jié)果，見表3.

表3 逐次分半辛普森算法計算結(jié)果Table.3 Results by using composite Simpson’s formula

由（9）式計算得

即身高中等（165cm，175cm）的大學生約為54.67%，不足60%．辛普森求積公式計算結(jié)果的分布曲線見圖3。

圖1 被積函數(shù)φ（x）的分布Fig.1 Distribution of the integrand

圖2 逐次分半梯形算法結(jié)果的分布Fig.2 Distribution of trapezoidal integration

圖3 逐次分半辛普森算法結(jié)果的分布Fig.3 Distribution of Simpson’s integration

圖4 三種算法的結(jié)果比較Fig.4 Comparison of three methods

將密度函數(shù)φ（x）的曲線、逐次分半梯形法計算結(jié)果的曲線、逐次分半辛普森法計算結(jié)果的曲線畫在同一坐標系下得比較圖，即圖4.由分布曲線可知，逐次分半梯形求積公式、逐次分半辛普森求積公計算結(jié)果基本符合正態(tài)分布，且分布曲線與標準分布曲線φ（x）擬合得很好.由表2和表3也可以看出，逐次分半辛普森公式收斂速度較快.

4.3 儲量計算問題

某地區(qū)為了估計某種礦物的儲量，考察隊來到該區(qū)進行勘測，得到如下數(shù)據(jù)（見表4）.請估計出此地區(qū)內(nèi)（1≤X≤4，1≤Y≤5）該礦物的儲量.

表4 地區(qū)勘探數(shù)據(jù)表Table.4 regional exploration data

圖5 數(shù)據(jù)點位置示意圖Fig.5 Locations of the data points

圖6 礦物的儲量分布Fig.6 Distribution of mineral reserves

利用公式（4）和（5），其中hx=hy=1，n=3，m=4，相應(yīng)的矩陣

經(jīng)計算，結(jié)果為256.320．統(tǒng)一單位后，得到該礦物的儲量約為2.56320×108（立方米），礦物儲量擬合圖見圖6.

由于給定的節(jié)點關(guān)于橫坐標x是奇數(shù)等分的分點，所以選用二重復合梯形求積公式計算較為合理．如果節(jié)點關(guān)于橫坐標x和縱坐標y都是偶數(shù)等分的分點，則也可以用二重復合辛普森求積公式計算.

5 結(jié)論

通過對上述三個實際問題的求解，不難發(fā)現(xiàn)，對于只給出有限個節(jié)點及其函數(shù)值的積分問題，針對不同的數(shù)據(jù)應(yīng)采用不同的分析、求解方法．一般情況下，逐次分半辛普森法的收斂速度是比較快的，甚至有時與龍貝格算法所需分半次數(shù)一樣.對于二重數(shù)值積分問題，如果被積函數(shù)的解析表達式未知，可利用復化梯形公式求積分的近似值.

［1］徐利治.數(shù)值積分法研究近況綜述［J］.數(shù)學進展，1982,11（2）：135-143.

［2］李慶揚等.《數(shù)值分析》第五版［M］.北京：清華大學出版社，2008.12：97-102

［3］王少英.數(shù)值積分若干方法的比較分析［J］.2012,28（6）：14-18.

［4］張智豐.《數(shù)學實驗》［M］.北京：科學出版社，2008.

［5］NENAD UJEVIC.Double integral inequalities and application in numerical integration［J］.Periodica Mathematica Hungarica,2004,49（1）：141-149.

［6］董志遠等.二重數(shù)值積分的計算方法的比較［J］.太原師范學院學報（自然科學版），2009,8（4）：20-23.

［7］AVRAM SIDI.Extension of a class of periodizing variablerans formationst for numerical integration［J］MATHEMATICS OF COMPUTATION,2005,75（253）：327-343.

［8］蘇岐芳等.《數(shù)值分析》［M］.北京：中國鐵道出版社，2007.

Iterative Methods of Numerical Integration and its Applications

HE Jun-jun，SU Qi-fang＊

（School of Mathematics and Information Engineering,Taizhou university,Linhai 317000,China）

We often encounter various types of integration problems in engineering and scientific computation.The integrand is too complex, even not the elementary function; or the analytical formula of the integrand is unknown.Only some nodes are given.There is a need to solve the approximation of integration.In this paper,we give two iterative algorithms of numerical integration and composite trapezoidal role for double integration, and solve some practical problems by using these iterative methods.

Numerical integration;algorithm;convergence speed；MATLAB program

10.13853/j.cnki.issn.1672-3708.2014.03.001

2014-05-05；

2014-05-20

簡介：蘇岐芳（1964- ），女，黑龍江綏化人，副教授，碩士，主要從事計算數(shù)學與應(yīng)用數(shù)學的教學與研究。