各種數(shù)值積分

1、各種數(shù)值積分第1頁，共85頁，2022年，5月20日，7點25分，星期二5.1 引言 利用牛頓萊布尼茲（NewtonLeibniz）公式 (5.1) 解決函數(shù) 在 上的積分問題在理論和應用上都有重大的意義。然而，在實際問題中，往往會遇到一些困難。有些形式上較簡單的函數(shù)，其原函數(shù) 不易求出或不能用初等函數(shù)表示成有限形式；有些被積函數(shù)的原函數(shù)過于復雜；而有些函數(shù)的函數(shù)值是由實驗、觀測等方法得出，并沒有給出具體的解析表達式。這些情形說明公式(5.1)在應用上是有局限性的，因此研究定積分的數(shù)值計算問題就顯得十分必要。 本章主要介紹一些常用的數(shù)值積分方法，包括梯形積分法、辛卜生積分法、變步長積分法、牛頓

2、柯特斯積分法、高斯積分法、龍貝格積分法以及高振蕩函數(shù)的積分法。第2頁，共85頁，2022年，5月20日，7點25分，星期二5.2 梯形積分法 5.2.1 梯形積分法的基本思想 梯形積分法的基本思想：在積分區(qū)間 上，根據(jù)給定的插值條件 和 ，構(gòu)造一個一次二項式 ，并以 的積分值近似地代替 。從幾何角度而言，是以梯形面積近似地代替曲邊梯形的面積。第3頁，共85頁，2022年，5月20日，7點25分，星期二5.2.2 梯形求積公式 依據(jù)梯形積分法的基本思想，將區(qū)間 分成 個 相等的小區(qū)間,則每個小區(qū)間的長度為 ,對每個小區(qū)間均實施如下的梯形求積： 將這些小梯形的求積值加起來，可以得到如下梯形求積公式

3、： 其中，第4頁，共85頁，2022年，5月20日，7點25分，星期二5.2.3 實現(xiàn)梯形積分法的基本步驟 (1) 輸入?yún)^(qū)間 的端點 值以及分割數(shù) ； (2) 將區(qū)間 等分成 個小區(qū)間，每一個小區(qū)間的 長度 ； (3) 計算每一個等分點的函數(shù)值 (4) 計算 (5) 輸出 的值； (6) 結(jié)束。 第5頁，共85頁，2022年，5月20日，7點25分，星期二圖 5.2 梯形積分法的N-S圖描述 第6頁，共85頁，2022年，5月20日，7點25分，星期二例5.1 使用梯形求積公式求下列定積分的值。第7頁，共85頁，2022年，5月20日，7點25分，星期二#define N 16 /* 等分數(shù)

4、*/float func(float x) float y; y=4.0/(1+x*x); return(y);void gedianzhi(float y,float a,float h) int i; for(i=0;i=N;i+) yi=func(a+i*h);float trapeze(float y,float h) float s; int i; s=(y0+yN)/2.0; for(i=1;iN;i+) s+=yi; return(s*h);第8頁，共85頁，2022年，5月20日，7點25分，星期二main() float a,b,h,s,fN+1; clrscr(); pri

5、ntf(input a,b=); scanf(%f,%f,&a,&b); h=(b-a)/(float)N; gedianzhi(f,a,h); s=trapeze(f,h); printf(s=%fn,s);程序運行結(jié)果：input a,b=0,1s=3.140942第9頁，共85頁，2022年，5月20日，7點25分，星期二 辛卜生積分法的基本思想：在積分區(qū)間 上，根據(jù)給定的插值條件 、 和 ，構(gòu)造個二次插值求積多項式 ，并以 的積分值近似地代替 。從幾何角度而言，是用過三點的拋物線面積近似地代替積分的曲邊面積。 5.3 辛卜生（Simpson）積分法 5.3.1 辛卜生積分法的基本思想第

6、10頁，共85頁，2022年，5月20日，7點25分，星期二5.3.2 辛卜生求積公式 依據(jù)辛卜生積分法的基本思想，將區(qū)間 分成 ( 必須是偶數(shù)）個相等的小區(qū)間，則每個小區(qū)間的長度為 ，在小區(qū)間 均實施如下的辛卜生求積： 將這些求積值加起來，可以得到如下辛卜生求積公式： 其中: 為寄數(shù)項的函數(shù) 值之和。為偶數(shù)項的函數(shù) 值之和。第11頁，共85頁，2022年，5月20日，7點25分，星期二5.3.3 實現(xiàn)辛卜生積分法的基本步驟(1) 輸入?yún)^(qū)間 的端點的 值以及分割數(shù) ；(2)將區(qū)間 等分成 個小區(qū)間，每一個小區(qū)間的長度 ；(3) 計算每一個等分點的函數(shù)值 ；(4) 計算: （計算奇數(shù)項的函數(shù)值之

7、和） （計算偶數(shù)項的函數(shù)值之和）(5) 計算 ；(6) 輸出 的值；(7) 結(jié)束。第12頁，共85頁，2022年，5月20日，7點25分，星期二圖 5.4 辛卜生積分法的N-S圖描述第13頁，共85頁，2022年，5月20日，7點25分，星期二例5.2 使用辛卜生求積公式求下列定積分的值。第14頁，共85頁，2022年，5月20日，7點25分，星期二#include #define N 16 /* 等分數(shù) */float func(float x) float y; y=4.0/(1+x*x); return(y);void gedianzhi(float y,float a,float h)

8、 int i; for(i=0;i=N;i+) yi=func(a+i*h);第15頁，共85頁，2022年，5月20日，7點25分，星期二float simpson(float y,float h) float s,s1,s2; int i; s1=y1; s2=0.0; for(i=2;i=eps) /* 判斷是否達到精度要求,若沒有達到,繼續(xù)循環(huán) */ s=0.0; for(i=0;i=n-1;i+) x=a+(i+0.5)*h; s=s+func(x); t2=(t1+h*s)/2.0; /* 計算 */ p=fabs(t1-t2); /* 計算精度 */ t1=t2; n=n+n;

9、h=h/2.0; return(t2);第22頁，共85頁，2022年，5月20日，7點25分，星期二void main() float a,b,t; clrscr(); printf(input a,b=); scanf(%f,%f,&a,&b); t=btrapeze(a,b); printf(t=%fn,t);程序運行結(jié)果：input a,b=0,1t=0.746824第23頁，共85頁，2022年，5月20日，7點25分，星期二 5.4.4 變步長辛卜生求積分法變步長辛卜生求積分法的基本過程:（1）利用梯形公式，將積分區(qū)間 一等分，（2）將其中的每一個求積小區(qū)間再二等分一次（3）根據(jù)上

10、面兩式 和 的余項 、 ，可以 推導出如下的變步長辛卜生求積公式。 進一步得到再二次等分一次后的變步長辛卜生求積公式為（4）若 ，二等分后的積分值 就是最后的結(jié)果；否則保存當前的變步長梯形積分值、等分數(shù)、積分值與步長，轉(zhuǎn)到第（2）步繼續(xù)做二等分處理。第24頁，共85頁，2022年，5月20日，7點25分，星期二5.4.5 實現(xiàn)變步長辛卜生積分法的基本步驟第25頁，共85頁，2022年，5月20日，7點25分，星期二第26頁，共85頁，2022年，5月20日，7點25分，星期二變步長梯形積分法的N-S圖描述第27頁，共85頁，2022年，5月20日，7點25分，星期二例5.4 使用變步長辛卜生求

11、積分法求下列定積分的值。#include #include #define eps 0.000001 /* 容許誤差 */float func(float x) float y; y=sqrt(1-x*x); return(y);第28頁，共85頁，2022年，5月20日，7點25分，星期二 float bsimpson(float a,float b) int i,n; float h,p,e,s; float t1,t2,s1,s2,x; n=1; h=b-a; t1=h*(func(a)+func(b)/2.0; s1=t1; /* 用代替 */ e=eps+1.0; while(e=e

12、ps) s=0.0; for(i=0;i=n-1;i+) x=a+(i+0.5)*h; s=s+func(x); t2=(t1+h*s)/2.0; /* 計算 */ s2=(4*t2-t1)/3.0; /* 計算 */ e=fabs(s2-s1); /* 計算精度 */ t1=t2; s1=s2; n=n+n; h=h/2.0; return(s2); 第29頁，共85頁，2022年，5月20日，7點25分，星期二void main() float a,b,s; clrscr(); printf(input a,b=); scanf(%f,%f,&a,&b); s=bsimpson(a,b);

13、 printf(s=%fn,s);程序運行結(jié)果：input a,b=0,1s=0.785398第30頁，共85頁，2022年，5月20日，7點25分，星期二5.5 牛頓柯特斯(NewtonCotes)積分法 5.5.1 牛頓柯特斯積分法的基本思想 牛頓柯特斯積分法的基本思想：用高次的插值求積多項式 去逼近被積函數(shù) ，以獲得高精度的積分值。 事實上，梯形積分是當 時的牛頓柯特斯積分，辛卜生積分是當 時的牛頓柯特斯積分，它們都是牛頓柯特斯積分的特例。第31頁，共85頁，2022年，5月20日，7點25分，星期二5.5.2 牛頓柯特斯求積公式 下面給出三到五階牛頓柯特斯求積公式。第32頁，共85頁，

14、2022年，5月20日，7點25分，星期二實現(xiàn)三階牛頓柯特斯求積公式的基本步驟如下：第33頁，共85頁，2022年，5月20日，7點25分，星期二三階牛頓柯特斯求積公式的N-S圖描述第34頁，共85頁，2022年，5月20日，7點25分，星期二例5.5 使用牛頓柯特斯求積公式求下列定積分的值。第35頁，共85頁，2022年，5月20日，7點25分，星期二#include #define N 5float func(float x) float y; y=4.0/(1+x*x); return(y);void gedianzhi(float y,float a,float b,int n) fl

15、oat h,s; int i; h=(b-a)/(float)n; for(i=0;i=n;i+) yi=func(a+i*h);第36頁，共85頁，2022年，5月20日，7點25分，星期二float nc3(float y,float a,float b,int n) float h,s,s0,s1; int i; h=(b-a)/(float)n; s0=s1=0.0; for(i=0;i=n-3;i=i+3) s0+=yi+yi+3; s1+=yi+1+yi+2; s=s0+3.0*s1; return(s*h*3.0/8.0);main() float a,b,h,s,fN+1; f

16、loat n3,n4,n5; printf(input a,b=); scanf(%f,%f,&a,&b); gedianzhi(f,a,b,3); n3=nc3(f,a,b,3); printf(n 3-nc 4-nc 5-ncn); printf(%f %f %fn,n3,n4,n5);第37頁，共85頁，2022年，5月20日，7點25分，星期二程序運行結(jié)果：input a,b=0,1 3-nc 4-nc 5-nc 3.138461 3.142118 3.141878第38頁，共85頁，2022年，5月20日，7點25分，星期二5.6 高斯積分法 5.6.1 高斯積分法的基本思想 前面介

17、紹的幾種數(shù)值積分方法，都是先尋找一個插值求積多項式 ，并以 近似代替函數(shù) 進行積分而獲得積分的近似值，即 (5.2) 表明，函數(shù) 在區(qū)間 上的積分可以用函數(shù) 在該區(qū)間上 的 個點的函數(shù)值的線性組合來近似代替。但是，由于插值求積公式是利用插值多項式的積分得到的，因此，如果被積函數(shù) 為次數(shù)不超過 的多項式，則利用插值求積公式計算得到的積分值是準確的。 在實際應用中，為了提高數(shù)值求積公式的精度，一般要求數(shù)值求積公式對于次數(shù)盡可能高的多項式能準確成立。由此提出了高斯積分法。即：如果插值求積公式(5.2)具有 次代數(shù)精度，那么稱該插值求積公式為高斯求積公式。其節(jié)點 稱為高斯點， 稱為高斯求積系數(shù)。 下面

18、具體介紹幾個常用的高斯求積公式。 第39頁，共85頁，2022年，5月20日，7點25分，星期二5.6.2 勒讓德高斯(Legendre-Gauss)求積公式 勒讓德高斯求積公式，特別適合于計算區(qū)間-1,1的積分，其求積公式可以表示為以下形式： 而對于一般區(qū)間 ，通過變換可以得到以下形式的求積公式： 高斯點 及高斯求積系數(shù) ，參見表5.1（ 表示階數(shù)）。第40頁，共85頁，2022年，5月20日，7點25分，星期二第41頁，共85頁，2022年，5月20日，7點25分，星期二第42頁，共85頁，2022年，5月20日，7點25分，星期二第43頁，共85頁，2022年，5月20日，7點25分，星

19、期二第44頁，共85頁，2022年，5月20日，7點25分，星期二第45頁，共85頁，2022年，5月20日，7點25分，星期二第46頁，共85頁，2022年，5月20日，7點25分，星期二第47頁，共85頁，2022年，5月20日，7點25分，星期二第48頁，共85頁，2022年，5月20日，7點25分，星期二例5.6 使用勒讓德高斯求積公式求下列定積分的值。第49頁，共85頁，2022年，5月20日，7點25分，星期二#include double func(double x) double y; y=x*x+sin(x); return(y);double legendre_gauss(

20、double a,double b,int m,int n) double h,hx,y,s,dx,x0; int i,k; static double x=-0.9061798459,-0.5384693101,0.0000000000, 0.5384693101,0.9061798459; static double w=0.2369268851,0.4786286705,0.5688888889, 0.4786286705,0.2369268851; dx=(b-a)/(double)m; hx=dx/2; s=0.0; for(k=0;km;k+) x0=a+(double)k*dx+

21、hx; for(i=0;in;i+) y=x0+xi*hx; s=s+wi*func(y); return(s*hx);第50頁，共85頁，2022年，5月20日，7點25分，星期二main() double a,b,s; int m,n; clrscr(); printf(input duandian a,b=); scanf(%f,%f,&a,&b); printf(input jieshu m=); scanf(%d,&m); printf(input fengeshu n=); scanf(%d,&n); s=legendre_gauss(a,b,m,n); printf(s=%lfn

22、,s);程序運行結(jié)果：input duandian a,b=2.5,8.4input jieshu m=10input fengeshu n=5s=192.077774第51頁，共85頁，2022年，5月20日，7點25分，星期二5.6.3 埃爾米特高斯(Hermite-Gauss)求積公式 埃爾米特高斯求積公式，特別適合于計算如下形式的積分： 其求積公式可以表示為以下形式： 高斯點 及高斯求積系數(shù) ，參見表5.2（ 表示階數(shù)）。第52頁，共85頁，2022年，5月20日，7點25分，星期二第53頁，共85頁，2022年，5月20日，7點25分，星期二第54頁，共85頁，2022年，5月20日

23、，7點25分，星期二埃爾米特高斯求積公式的N-S圖描述 第55頁，共85頁，2022年，5月20日，7點25分，星期二例5.7 使用埃爾米特高斯求積公式求下列定積分的值。第56頁，共85頁，2022年，5月20日，7點25分，星期二#include double func(double x) double y; y=x*x; return(y);double hermite_gauss(int n) double s; int i,k; static double x=-2.02018287046,-0.95857246461,0.00000000000, 0.95857246461,2.02

24、018287046; static double w=0.01995324206,0.39361932315,0.94530872048, 0.39361932315,0.01995324206; s=0.0; for(i=0;in;i+) s=s+wi*func(xi); return(s);第57頁，共85頁，2022年，5月20日，7點25分，星期二main() double s; int n; clrscr(); printf(input n=); scanf(%d,&n); s=hermite_gauss(n); printf(s=%lfn,s);程序運行結(jié)果：input n=5s=

25、 0.886227第58頁，共85頁，2022年，5月20日，7點25分，星期二5.6.4 拉蓋爾高斯(Laguerre-Gauss)求積公式 拉蓋爾高斯求積公式，特別適合于計算如下形式的積分： 其求積公式可以表示為以下形式： 高斯點 及高斯求積系數(shù) ，參見表5.3（ 表示階數(shù)）。第59頁，共85頁，2022年，5月20日，7點25分，星期二第60頁，共85頁，2022年，5月20日，7點25分，星期二第61頁，共85頁，2022年，5月20日，7點25分，星期二拉蓋爾高斯求積公式的N-S圖描述第62頁，共85頁，2022年，5月20日，7點25分，星期二例5.8 使用拉蓋爾-高斯求積公式求下

26、列定積分的值。第63頁，共85頁，2022年，5月20日，7點25分，星期二#include double func(double x) double y; y=x; return(y);double hermite_gauss(int n) double s; int i,k; static double x=0.26356031972,1.41340305911,3.59642577104, 7.08581000586,12.64080084428; static double w=0.52175561058,0.39866681108,0.07594244968, 0.003611758

27、68,0.00002336997; s=0.0; for(i=0;in;i+) s=s+wi*func(xi); return(s); 第64頁，共85頁，2022年，5月20日，7點25分，星期二main() double s; int n; clrscr(); printf(input n=); scanf(%d,&n); s=hermite_gauss(5); printf(s=%lfn,s);程序運行結(jié)果：input n=5s= 1.000000第65頁，共85頁，2022年，5月20日，7點25分，星期二5.7 龍貝格（Romberg）積分法 5.7.1 龍貝格積分法的基本思想 前面

28、講述的各種求積方法是插值求積的思想，而龍貝格積分法的基本思想是，使用一個諸如梯形求積法等代數(shù)精度較低的求積公式，相繼以步長 和 求得定積分的兩個近似結(jié)果，然后再做它們適當?shù)木€性組合，就可以得到一個代數(shù)精度更高的公式。 第66頁，共85頁，2022年，5月20日，7點25分，星期二5.7.2 實現(xiàn)龍貝格積分法的基本步驟（1）輸入?yún)^(qū)間 的端點 的值,最大迭代次數(shù) 以及容許誤差；（2）計算區(qū)間 的長度 ；（3）用梯形積分法計算積分近似值 ；（4）對 計算 對 計算 ， 如果 ，則退出循環(huán)。（5）如果 ，則繼續(xù)；否則輸出無解信息，轉(zhuǎn)（7）；（6）輸出 的值；（7）結(jié)束。第67頁，共85頁，2022年，

29、5月20日，7點25分，星期二 表5.1龍貝格求積算法元素進行運算的順序?qū)崿F(xiàn)龍貝格積分法的NS圖，如圖5.11所示。 第68頁，共85頁，2022年，5月20日，7點25分，星期二例5.9 使用龍貝格求積公式求下列定積分的值。第69頁，共85頁，2022年，5月20日，7點25分，星期二#include #include #define DFS_N 20 /* 等分數(shù) */#define MAX_N 10 /* 最大循環(huán)次數(shù) */#define eps 0.00001 /* 容許誤差 */double func(double x) double y; y=4.0/(1+x*x); return

30、(y);double sum(double aa,double bb,long int n) double h,s; int i; h=(bb-aa)/n; s=0.0; for(i=1;in;i+) s+=func(aa+i*h); s=s+(func(aa)+func(bb)/2.0; return(h*s);第70頁，共85頁，2022年，5月20日，7點25分，星期二void romberg(double a,double b) double s,tMAX_N+12; int i,flag=0; long int n=DFS_N,m; t01=sum(a,b,n); n*=2; for

31、(m=1;mMAX_N;m+) for(i=0;im;i+) ti0=ti1; t01=sum(a,b,n); n*=2; for(i=1;i=m;i+) ti1=ti-11+(ti-11-ti-10)/(pow(2,2*m)-1); if(fabs(tm1-tm-11)eps) printf(t%ld0=%lfn,m,tm1); flag=1; break; if(flag=0) printf(Return no solovedn);第71頁，共85頁，2022年，5月20日，7點25分，星期二main() double a,b; clrscr(); printf(input a,b=);

32、scanf(%lf,%lf,&a,&b); romberg(a,b);程序運行結(jié)果：input a,b=0,1t20=3.141570第72頁，共85頁，2022年，5月20日，7點25分，星期二5.8 高振蕩函數(shù)的積分法 5.8.1 高振蕩函數(shù)的積分法的基本思想 在工程實際問題中，經(jīng)常會遇到如下形如 的積分，當m充分大時為高振蕩函數(shù)的積分。對于高振蕩函數(shù)的積分，如果采用插值求積法進行積分，則在建立被積函數(shù) 或 的插值多項式 時，為了使 能夠很好地逼近它們，就要求 也要振蕩得厲害，即要求插值多項式 的次數(shù)足夠高。但是，高次插值實際的逼近性質(zhì)很不好，實用價值不大。即使采用分段低次插值，效果也不會

33、很理想。 因此，引進計算高振蕩函數(shù)的積分的重要方法分部積分法。 第73頁，共85頁，2022年，5月20日，7點25分，星期二分部積分法的基本思想是， 令： 其中： 則有： 反復利用分部積分法可以得到 分離出實部和虛部后就得到以下分部積分公式。第74頁，共85頁，2022年，5月20日，7點25分，星期二5.8.2 分部積分公式當積分區(qū)間為 時，則變?yōu)?第75頁，共85頁，2022年，5月20日，7點25分，星期二實現(xiàn)高振蕩函數(shù)積分的NS圖第76頁，共85頁，2022年，5月20日，7點25分，星期二例5.10 用分部積分法計算下列高振蕩積分的值。 其中 ， ， 。取 ， ， ， 則有 第77

34、頁，共85頁，2022年，5月20日，7點25分，星期二#include void part(double a,double b,int m,int n,double fa,double fb,double s) int mm,i,j; double sma,smb,cma,cmb; double sa4,sb4,ca4,cb4; sma=sin(m*a); smb=sin(m*b); cma=cos(m*a); cmb=cos(m*b); sa0=sma; sa1=cma; sa2=-sma; sa3=-cma; sb0=smb; sb1=cmb; sb2=-smb; sb3=-cmb; ca0=cma; ca1=-sma; ca2=-cma; ca3=sma; cb0=cmb; cb1=-smb; cb2=-cmb; cb3=smb; s0=0.0; s1=0.0;