《数值分析》课程教学资源（PPT讲稿）第八章 数值积分（2/2）

§2复合求积/ Composite Quadrature 高次插值有 Runge现象,故采用分段低次插值 →分段低次合成的 Newton-Cotes复合求积公式 >复合梯形公式:h b-a xk=a+kh(k=0,…,n) 在每个x1x上用梯形公式 88 (sc:3pmP如野yx),k= k-1 e oscllator ∫/么gP=7岛+2)+/b k=1 ∑∫(5k R=∑ h f"(5k=-(b-a /中值定理* 12 h (b-a)f"(95),5∈(a,b)

§2 复合求积 /* Composite Quadrature */ Haven’t we had enough formulae? What’s up now? Oh come on, you don’t seriously consider h=(b−a)/2 acceptable, do you? Why can’t you simply refine the partition if you have to be so picky? Don’t you forget the oscillatory nature of high￾degree polynomials! Uh-oh 高次插值有Runge 现象，故采用分段低次插值  分段低次合成的 Newton-Cotes 复合求积公式。 ➢ 复合梯形公式： , x a k h (k 0, ... ,n) n b a h k = + = − = 在每个 [xk−1 , xk ] 上用梯形公式： f x f x k n x x f x dx k k x x k k k k [ ( ) ( )], 1, ... , 2 ( ) 1 1 1 + = −  − −  −       = +  + − = 1 1 ( ) 2 ( ) ( ) 2 n k k f a f x f b h  =  − + b a n k k k f x f x h f x dx 1 1 [ ( ) ( )] 2 ( ) = Tn ( ) ( ), ( , ) 12 ( ) ( ) 12 ( )] 12 [ ] [ 2 1 2 1 3 b a f a b h n f b a h f h R f n k n k k k = − −    = −  = − −   = =     /*中值定理*/

8 2 Composite Quadrature >复化 Simpson公式:h=b-a,x4=a+kh(k=0,,m f(x)daf(xx)+4f(x+f(akI k k+1 ∫(x)k=/(+2(x+2(xm)+(=n R[/b-a/h)4 5) 180(2 注:为方便编程,可采用另一记法:令n’=2n为偶数 这时n=b=a=B,x1=a+kM,有 Sn=,Uf(a)+4∑f(x)+2∑f(xk)+f() odd k even k

§2 Composite Quadrature ➢ 复化 Simpson 公式： , x a k h (k 0, ... ,n) n b a h k = + = − = [ ( ) 4 ( ) ( )] 6 ( ) 1 2 1 1   + + + + + k k k x x f x f x f x h f x dx k k xk 2 1 k+ x xk+1 4 4 4 4 4 [ ( ) 4 ( ) 2 ( ) ( )] 6 ( ) 1 0 1 0 1 2   1  − = − =  + + + + + n k n k k k b a f a f x f x f b h f x dx = Sn ( ) 180 2 [ ] (4) 4 f  b a h R f      −  = − 注：为方便编程，可采用另一记法：令 n’ = 2n 为偶数， 这时 h x a k h ，有 n b a h k = = +   −  = , 2 [ ( ) 4 ( ) 2 ( ) ( )] 3 +  +  +  = odd k even k n k k f a f x f x f b h S

8 2 Composite Quadrature >收敛速度与误差估计: 定义若一个积分公式的误差满足hm刚c<a且C≠0, 则称该公式是p阶收敛的。 TnNOCh', s,o(h), cn och) 例:计算z 1+x 运算量基 本相同 解:16/(0)+2∑f(x)+/1)其中x= 3.138988494 S10+2()+2+0其中x 3.141592502

§2 Composite Quadrature ➢ 收敛速度与误差估计： 定义 若一个积分公式的误差满足 且C  0， 则称该公式是 p 阶收敛的。 =   → C h R f p h [ ] lim 0 ( ) , ( ) , ( ) 2 4 6 Tn ~ O h Sn ~ O h Cn ~ O h 例：计算 dx  + x = 1 0 1 4  2 解：       = +  + = (0) 2 ( ) (1) 16 1 7 1 8 T f f x f k k 8 k x 其中 k = = 3.138988494       = (0) + 4 ( ) + 2 ( ) + (1) 24 1 odd even 4 S f f x f x f k k 8 k x 其中 k = = 3.141592502 运算量基 本相同

82 Composite Quadrature Q:给定精度,如何取n? 例如:要求|Ⅰ-Tn|<E,如何判断n=? HW: p74 R∫]= (b-a)f"(5 ∑Ir5)川#2 h ∫"(x)x 12 Ja 上例中若要求-r<10,则|R/=h|r()-r(0)=b2<10 →h<0.00244949即:取n=409 通常采取将区间不断对分的方法,即取n=k 上例中2k≥409→k=9时,T5 可用来判断迭代 是否停止 注意到区间再次对分时R2f≈ 1592504R, L

§2 Composite Quadrature Q: 给定精度 ，如何取 n ? 例如：要求 | I −Tn |   ，如何判断 n = ？ ( ) ( ) 12 [ ] 2 b a f  h R f = − −  ? = = −   n k f k h h 1 2 [ ( ) ] 12  [ ( ) ( )] 12 ( ) 12 2 2 f b f a h f x dx h b a  −  = −  −   上例中若要求 ，则 6 | | 10− I −Tn  6 2 2 10 6 | (1) (0)| 12 | [ ]| −   −  =  h f f h R f n h 0.00244949 即：取 n = 409 通常采取将区间不断对分的方法，即取 n = 2k 上例中2 k  409  k = 9 时，T512 = 3.14159202 注意到区间再次对分时 [ ] 4 1 2 [ ( ) ( )] 12 1 [ ] 2 2 R f h R f f b f a n   n       −  −  4 2 1  − − n n I T I T ( ) 3 1 T2n T2n Tn I −  − 可用来判断迭代 是否停止。 HW: p.174 #2

8 2 Composite Quadrature Lab 13. Composite Trapezoidal rule Use the Composite Trapezoidal rule with a given n>0 to approximate a given integral f(x)dx You are supposed to write a function double Ctr( int n, double a, double b, double("Do) to approximate the integral from a to b of the function f using the trapezoidal rule on n equal-length subintervals. Input There is no input file. Instead, you must hand in your function in a*.h file. The rule of naming the *. h file is the same as that of naming the *. c or cpp files. Output For each test case, you are supposed to return the approximation of the integral

§2 Composite Quadrature Lab 13. Composite Trapezoidal Rule Use the Composite Trapezoidal rule with a given n > 0 to approximate a given integral . You are supposed to write a function double CTR ( int n, double a, double b, double (*f)( ) ) to approximate the integral from a to b of the function f using the trapezoidal rule on n equal-length subintervals. Input There is no input file. Instead, you must hand in your function in a *.h file. The rule of naming the *.h file is the same as that of naming the *.c or *.cpp files. Output For each test case, you are supposed to return the approximation of the integral.  b a f ( x) dx

8 2 Composite Quadrature Sample Judge Program #include fclose(outHile); #include # include"9811500113.h" double fl( double x) i return(1.0/(1.0+sin(x) *sin(x)));) Sample output double f2( double x) 0.809093 i return(x*log(x));) 0.639900 void main( f FILE* outhile= fopen ("out. txt",w; Int n: double a, b; 0;b=1.0;n=10: fprintf(outfile, %lf\n",CtR(, a, b, fD)) a=1.0;b=20;n=4 fprintf(outfile, %lfn", CTR(n, a, b, f2))

§2 Composite Quadrature Sample Judge Program #include #include #include "98115001_13.h" double f1 ( double x ) { return (1.0/(1.0+sin(x)*sin(x))); } double f2 ( double x ) { return (x*log(x)); } void main( ) { FILE *outfile = fopen("out.txt", "w"); int n; double a, b; a = 0.0; b = 1.0; n = 10; fprintf(outfile, "%lf\n", CTR(n, a, b, f1)); a = 1.0; b = 2.0; n = 4; fprintf(outfile, "%lf\n", CTR(n, a, b, f2)); fclose(outfile); } Sample Output 0.809093 0.639900

§3龙贝格积分 Romberg Integration 例:计算n=J。1 已知对于=10须将区间对分9次,得到T512=314159202 考察,a!由ls41nh=4rn-1r,来计算/效果是否好些? T-T4=3.141592502=S4 Romberg E列 一般有: IT-T =S 42S. =S C-C =C : R 42-1 <E? Romberg S1=T0<<y 算法:④T=r( GS2=T1、N⑥C1=T2 (0) ⑧S4=T ⑩R1=T2

§3 龙贝格积分 /* Romberg Integration */ 例：计算 dx  + x = 1 0 1 4  2 已知对于 = 10−6须将区间对分 9 次，得到 T512 = 3.14159202 由 n n T n Tn 来计算 I 效果是否好些？ T T I 3 1 3 4 4 1 4 2 2 = − − − 考察  4 1 2  − − n n I T I T 8 4 3 1 3 4 T − T = 3.141592502 = S4 一般有： n n n S T T = − − 4 1 4 2 n n n C S S = − − 4 1 4 2 2 2 n n n R C C = − − 4 1 4 3 2 3 Romberg 序列 ➢ Romberg 算法： <  ? <  ? <  ? … … … … … …  T1 = (0) T0  T8 = (3) T0  T4 = (2) T0  T2 = (1) T0  S1 = (0) T1  R1 = (0) T3  S2 = (1) T1  C1 = (0) T2  C2 = (1)  S4 = T2 (2) T1

8 3 Romberg integration Lab 14. romberg Integration Use Romberg integration to approximate a given integral/6 f(r)de You are supposed to write a function double romberg(double a, double b, double(fo, double eps, int*k) to approximate the integral from a to b of the function f using using Romberg integration. Output when T o)-To)< eps Input There is no input fle. Instead, you must hand in your function in a"h file. The rule of naming the *h file is the same as that of naming the *c or cpp files. Output For each test case, you are supposed to return the approximation of the integral and save the totaliteration number in k The maximum iteration number is set to be MAx N. If the romberg sequence is not convergent after maX N iterations, return the current Maxn, and set k to be -1

§3 Romberg Integration Lab 14. Romberg Integration Use Romberg integration to approximate a given integral . You are supposed to write a function double Romberg (double a, double b, double (*f)( ) , double eps, int *k) to approximate the integral from a to b of the function f using using Romberg integration. Output when eps. Input There is no input file. Instead, you must hand in your function in a *.h file. The rule of naming the *.h file is the same as that of naming the *.c or *.cpp files. Output For each test case, you are supposed to return the approximation of the integral and save the total iterationnumber in k. The maximum iteration number is set to be MAX_N. If the Romberg sequence is not convergent after MAX_N iterations, return the current , and set k to be −1.  b a f ( x) dx − −  (0) 1 (0) Tk Tk (0) TMAX_N

8 3 Romberg integration Sample Judge Program #include fprintfoutfile, %If\n", R); #include #define maxn 1000 fclose( outfile); # nclude"9811500113h" double fo(double x i return(1.0/x); Sample output void main( i FILE *outfile=fopen ( out.txt","W); 41.098631 int ks double a, b, eps, r; a=1.0;b=3.0;eps=0.01; R=Romberg(a, b, f0, eps, &k); if(k==-1) fprintf(outfile, Maximum number of iterations exceeded. n); else fprintf(outfile, %d",k);

§3 Romberg Integration Sample Judge Program #include #include #define MAX_N 1000 #include "98115001_13.h" double f0(double x) { return (1.0/x); } void main( ) { FILE *outfile = fopen("out.txt", "w"); int k; double a, b, eps, R; a = 1.0; b = 3.0; eps = 0.01; R = Romberg (a, b, f0, eps, &k); if (k == −1) fprintf(outfile, "Maximum number of iterations exceeded.\n"); else fprintf(outfile, "%d ", k); fprintf(outfile, "%lf\n", R); fclose(outfile); } Sample Output 4 1.098631

8 3 Romberg Integration >理查德森外推法/ Richardson's extrapolation G利用低阶公式产生高精度的结果。(a1与h无关 设对于某一h≠0,有公式T0)近似计算某不知值l。由 Taylor展开得到:r0(h)-I=a1h+a2l2+a3h3+ 现将h对分,得:(当)-1=a1(+a2()+a3(分)+ Q:如何将公式精度由O(h)提高到O(h2)? 2T0()-T0(h) 3 a, h3 2-1 2 2 p175 2T0()-Th #4 即:T(h) 1+B,h+B,h+ 2-1 2T()-T(h 22-1 1+y1h+y2h+… →Tn(h)= 2"Tm1(2)-Tm-1(h =Ⅰ+b,hm+1+S,hm+2+… 2m-1

§3 Romberg Integration ➢ 理查德森外推法 /* Richardson’s extrapolation */ 利用低阶公式产生高精度的结果。 设对于某一 h  0，有公式 T0 (h) 近似计算某一未知值 I。由 Taylor展开得到： T0 (h) − I = 1 h + 2 h 2 + 3 h 3 + … i 与 h 无关 现将 h 对分，得： ( ) () ( ) ( ) ... 3 3 2 2 0 2 1 2 2 2 T h − I =  h + h + h + Q：如何将公式精度由 O(h) 提高到 O(h 2 ) ？ ... 4 3 2 1 2 1 2 ( ) ( ) 3 3 2 2 0 2 0 − = − − − − − I h h T T h h   即： ... 2 1 2 ( ) ( ) ( ) 3 2 2 1 0 2 0 1 = + + + − − = I h h T T h T h h   ( ) ... 4 2 3 T2 h = = I + 1 h + h + 2 1 2 ( ) ( ) 2 1 2 1 2 − T −T h h ( ) ... 2 2 1 = = + 1 + + m+ m+ Tm h I  h  h 2 1 2 ( ) ( ) 1 2 1 − − − − m m h m m T T h HW: p.175 #4