§ 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 复合 求积公式。
复合梯形公式,),.,,,0(,nkhkax
n
abh
k
在每个 上用梯形公式：],[ 1 kk xx?
nkxfxfxxdxxf kkxx kkk
k
,...,1,)]()([2)( 11
1
1
1
)()(2)(2
n
k
k bfxfaf
h
b
a
n
k
kk xfxf
hdxxf
1
1 )]()([2)(
= Tn
),(),()(
12
)(
)(
12
)](
12
[][
2
1
2
1
3
bafab
h
n
f
ab
h
f
h
fR
n
k
kn
k
k
/*中值定理 */
§ 2 Composite Quadrature
复化 Simpson 公式：
),.,,,0(,nkhkaxn abh k
)]()(4)([6)( 1
21
1
kkk
x
x xfxfxf
hdxxfk
k
kx
21?kx
1?kx
4
4
4
4
4
])()(2)(4)([6)( 1
0
1
0
121
n
k
n
k
kk
b
a bfxfxfaf
hdxxf = Sn
)(21 8 0][ )4(
4
fhabfR
注,为方便编程，可采用另一记法：令 n’ = 2n 为偶数，
这时,有
hkaxhn abh k,2
])()(2)(4)([3
ko d d ke v e n
kkn bfxfxfaf
hS
§ 2 Composite Quadrature
收敛速度与误差估计：
定义 若一个积分公式的误差满足 且 C? 0，
则 称该公式是 p 阶收敛 的 。
Ch fR ph ][lim 0
)(,)(,)( 642 hOChOShOT nnn ~~~
例,计算 dx
x
1
0 1
4
2?
解：
)1()(2)0(161 7
1
8 fxffT
k
k8
kx
k?
其中
= 3.138988494
)1()(2)(4)0(241
o d d e v e n
4 fxfxffS kk8kxk?
其中
= 3.141592502
运算量基本相同
§ 2 Composite Quadrature
Q,给定精度?，如何取 n?
例如：要求,如何判断 n =？ || nTI
)()(12][ 2?fabhfR
n
k
k hf
h
1
2
])([12?
)]()([12)(12 22 afbfhdxxfh ba
上例中若要求,则610|| nTI 622 10
6|)0()1(|12|][|
hffhfR
n
0 0 2 4 4 9 4 9.0?h 即：取 n = 409
通常采取将区间 不断对分 的方法，即取 n = 2k
上例中 2k? 409? k = 9 时,T512 = 3.14159202
注意到区间再次对分时 ][
4
1
2)]()([12
1][ 2
2 fR
hafbffR
nn
4
12?
n
n
TI
TI )(
3
1
22 nnn TTTI
可用来判断迭代是否停止。
HW,
p.174
#2
§ 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.
ba dxxf )(
§ 2 Composite Quadrature
Sample Judge Program
#include <stdio.h>
#include <math.h>
#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 */
例,计算 dx
x
1
0 1
4
2?
已知对于?= 10?6须将区间对分 9 次，得到 T512 = 3.14159202
由 来计算 I 效果是否好些？
nnnn TT
TTI
3
1
3
4
14
4
22
考察
4
12?
n
n
TI
TI
48 3
1
3
4 TT? = 3.141592502 = S4
一般有：
nnn S
TT?
14
4 2
nnn C
SS?
14
4
2
2
2
nnn R
CC?
14
4
3
2
3
Romberg 序列
Romberg
算法：
<
<
<
… … … … … …
T1 = )0(0T
T8 = )3(0T
T4 = )2(0T
T2 = )1(0T? S1 = )0(1T
R1 = )0(3T
S2 = )1(1T? C1 = )0(2T
C2 = )1(2T? S4 = )2(1T
§ 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 iteration 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
,and set k to be?1.
ba dxxf )(
)0( 1)0( kk TT
)0(MAX_NT
§ 3 Romberg Integration
Sample Judge Program
#include <stdio.h>
#include <math.h>
#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
§ 3 Romberg Integration
理查德森 外推法 /* Richardson’s extrapolation */
利用 低 阶公式产生 高 精度的结果。
设对于某一 h? 0，有公式 T0(h) 近似计算某一未知值 I。由
Taylor展开得到,T0(h)? I =?1 h +?2 h2 +?3 h3 + …
i 与 h 无关现将 h 对分，得,() () (),..)( 3232222120 hhhh IT
Q,如何将公式精度由 O(h) 提高到 O(h2)？
...432112 )()(2 3322020 hhIhTT
h
即：
...12 )()(2)( 32210201 hhIhTThT h
...)( 42312 hhIhT
12
)()(2
2
121
2
hTT h
...)( 2211 mmm hhIhT
12
)()(2 121
m
m
h
m
m hTT
HW,
p.175
#4
§ 4 高斯型 积分 /* Gaussian Quadrature */
ba nk kk xfAdxxfx 0 )()()(?
构造具有 2n+1次代数精度的求积公式将节点 x0 … xn 以及系数 A0 … An 都作为待定系数。
令 f (x) = 1,x,x2,…,x2n+1 代入可求解，得到的公式具有 2n+1 次代数精度。这样的节点称为 Gauss 点,
公式称为 Gauss 型求积公式 。
例,求 的 2 点 Gauss 公式。dxxfx )(1
0?
解,设,应有 3 次代数精度。1
0 1100 )()()( xfAxfAdxxf
x
代入 f (x) = 1,x,x2,x3
3
11
3
009
2
2
11
2
007
2
11005
2
103
2
xAxA
xAxA
xAxA
AA
2776.0
3891.0
2899.0
8212.0
1
0
1
0
A
A
x
x
不是线性方程组，
不易求解。
§ 4 Gaussian Quadrature
证明：,?” x0 … xn 为 Gauss 点,则公式至少有 2n+1 次代数精度。
ba n
k
kk xfAdxxfx
0
)()()(?
对任意次数 不大于 n的多项式 Pm(x),Pm(x) w(x)的次数不大于 2n+1,则代入公式应 精确成立,
n
k
kkmk
b
a m xwxPAdxxwxPx 0 )()()()()(?
0
= 0?
“?” 要证明 x0 … xn 为 Gauss 点，即要证公式对任意次数 不大于 2n+1的多项式 Pm(x) 精确成立，即证明：
n
k
kmk
b
a m xPAdxxPx 0 )()()(?
设 )()()()( xrxqxwxP
m
bababa m dxxrxdxxqxwxdxxPx )()()()()()()(
0
n
k
kk xrA
0
)(
n
k
kmk xPA
0
)(?
x0… xn 为 Gauss 点? 与任意次数不大于 n的多项式 P(x) （ 带权 ） 正交 。
n
k
kxxxw
0
)()(
定理求 Gauss 点? 求 w(x)
§ 4 Gaussian Quadrature
正交多项式族 {?0,?1,…,?n,… }有性质：任意次数不大于 n 的多项式 P(x) 必与?n+1正交。
若取 w(x) 为其中的?n+1，则?n+1的根 就是 Gauss 点。
再解上例,10 1100 )()()( xfAxfAdxxfx
Step 1,构造正交多项式?2
设 cbxxxaxxx 2210 )(,)(,1)(
53a0)(10 dxaxx0),( 10
1
021
1
0
2
10
0))(53(0),(
0)(0),(
dxcbxxxx
dxcbxxx
21
5
9
10
c
b
即：
21
5
9
10)( 2
2 xxx?
§ 4 Gaussian Quadrature
Step 2,求?2 = 0 的 2 个根，即为 Gauss 点 x0,x1
2
21/20)9/10(9/10 2
1;0
x
Step 3,代入 f (x) = 1,x 以求解 A0,A1
解 线性 方程组，
简单。
结果与前一方法相同,2776.0,3891.0,2899.0,8212.0
1010 AAxx
利用此公式计算 的值? 1
0 dxex
x
2555.1 1
0 dxex
x 2899.08212.010 2 7 7 6.03 8 9 1.010 eeeAeA xx
注,构造正交多项式也可以利用 L-S 拟合中介绍过的递推式进行。
§ 4 Gaussian Quadrature
特殊正交多项式族：
① Legendre 多项式族,1)(?x?定义在 [?1,1]上，
k
k
k
kk xdx
d
kxP )1(!2
1)( 2 满足：
lk
lkPP
k
lk
12
2
0),(
xPP 10,1由 有递推 11 )12()1( kkk kPxPkPk
以 Pn+1 的根为节点的求积公式称为 Gauss-Legendre 公式 。
② Chebyshev 多项式族,21 1)( xx定义在 [?1,1]上，
) a rcc o s( co s)( xkxT k
Tn+1 的根为
22 12c o s nkx kk = 0,…,n
以此为节点构造公式
1
1 02 )()(1
1 n
k
kk xfAdxxfx
称为 Gauss-Chebyshev 公式 。
注意到积分端点?1 可能是积分的 奇点,用普通 Newton-Cotes公式在端点会出问题。而 Gauss公式可能避免此问题的发生。
其它公式见教材 p.189
§ 4 Gaussian Quadrature
Gauss 公式的余项：
ba nk kk xfAdxxffR 0 )()(][ /* 设 P为 f的过 x0 … xn的插值多项式 */
ba nk kk xPAdxxf 0 )()(
/*只要 P 的阶数不大于 2n+1，则下一步等式成立 */
dxxPxfdxxPdxxf baba ba )]()([)()(
插值多项式 的余项Q,什么样的 插值多项式 在 x0 … xn 上有 2n+1 阶？
A,Hermite 多项式！ 满足 )()(),()(
kkkk xfxHxfxH
ba dxxHxffR )]()([][
),(,)(
)!22(
)(
)(
)!22(
)(
2
)12(
2
)12(
badxxw
n
f
dxxw
n
f
b
a
n
b
a
x
n
HW,
p.190 #6
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 复合 求积公式。
复合梯形公式,),.,,,0(,nkhkax
n
abh
k
在每个 上用梯形公式：],[ 1 kk xx?
nkxfxfxxdxxf kkxx kkk
k
,...,1,)]()([2)( 11
1
1
1
)()(2)(2
n
k
k bfxfaf
h
b
a
n
k
kk xfxf
hdxxf
1
1 )]()([2)(
= Tn
),(),()(
12
)(
)(
12
)](
12
[][
2
1
2
1
3
bafab
h
n
f
ab
h
f
h
fR
n
k
kn
k
k
/*中值定理 */
§ 2 Composite Quadrature
复化 Simpson 公式：
),.,,,0(,nkhkaxn abh k
)]()(4)([6)( 1
21
1
kkk
x
x xfxfxf
hdxxfk
k
kx
21?kx
1?kx
4
4
4
4
4
])()(2)(4)([6)( 1
0
1
0
121
n
k
n
k
kk
b
a bfxfxfaf
hdxxf = Sn
)(21 8 0][ )4(
4
fhabfR
注,为方便编程，可采用另一记法：令 n’ = 2n 为偶数，
这时,有
hkaxhn abh k,2
])()(2)(4)([3
ko d d ke v e n
kkn bfxfxfaf
hS
§ 2 Composite Quadrature
收敛速度与误差估计：
定义 若一个积分公式的误差满足 且 C? 0，
则 称该公式是 p 阶收敛 的 。
Ch fR ph ][lim 0
)(,)(,)( 642 hOChOShOT nnn ~~~
例,计算 dx
x
1
0 1
4
2?
解：
)1()(2)0(161 7
1
8 fxffT
k
k8
kx
k?
其中
= 3.138988494
)1()(2)(4)0(241
o d d e v e n
4 fxfxffS kk8kxk?
其中
= 3.141592502
运算量基本相同
§ 2 Composite Quadrature
Q,给定精度?，如何取 n?
例如：要求,如何判断 n =？ || nTI
)()(12][ 2?fabhfR
n
k
k hf
h
1
2
])([12?
)]()([12)(12 22 afbfhdxxfh ba
上例中若要求,则610|| nTI 622 10
6|)0()1(|12|][|
hffhfR
n
0 0 2 4 4 9 4 9.0?h 即：取 n = 409
通常采取将区间 不断对分 的方法，即取 n = 2k
上例中 2k? 409? k = 9 时,T512 = 3.14159202
注意到区间再次对分时 ][
4
1
2)]()([12
1][ 2
2 fR
hafbffR
nn
4
12?
n
n
TI
TI )(
3
1
22 nnn TTTI
可用来判断迭代是否停止。
HW,
p.174
#2
§ 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.
ba dxxf )(
§ 2 Composite Quadrature
Sample Judge Program
#include <stdio.h>
#include <math.h>
#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 */
例,计算 dx
x
1
0 1
4
2?
已知对于?= 10?6须将区间对分 9 次，得到 T512 = 3.14159202
由 来计算 I 效果是否好些？
nnnn TT
TTI
3
1
3
4
14
4
22
考察
4
12?
n
n
TI
TI
48 3
1
3
4 TT? = 3.141592502 = S4
一般有：
nnn S
TT?
14
4 2
nnn C
SS?
14
4
2
2
2
nnn R
CC?
14
4
3
2
3
Romberg 序列
Romberg
算法：
<
<
<
… … … … … …
T1 = )0(0T
T8 = )3(0T
T4 = )2(0T
T2 = )1(0T? S1 = )0(1T
R1 = )0(3T
S2 = )1(1T? C1 = )0(2T
C2 = )1(2T? S4 = )2(1T
§ 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 iteration 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
,and set k to be?1.
ba dxxf )(
)0( 1)0( kk TT
)0(MAX_NT
§ 3 Romberg Integration
Sample Judge Program
#include <stdio.h>
#include <math.h>
#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
§ 3 Romberg Integration
理查德森 外推法 /* Richardson’s extrapolation */
利用 低 阶公式产生 高 精度的结果。
设对于某一 h? 0，有公式 T0(h) 近似计算某一未知值 I。由
Taylor展开得到,T0(h)? I =?1 h +?2 h2 +?3 h3 + …
i 与 h 无关现将 h 对分，得,() () (),..)( 3232222120 hhhh IT
Q,如何将公式精度由 O(h) 提高到 O(h2)？
...432112 )()(2 3322020 hhIhTT
h
即：
...12 )()(2)( 32210201 hhIhTThT h
...)( 42312 hhIhT
12
)()(2
2
121
2
hTT h
...)( 2211 mmm hhIhT
12
)()(2 121
m
m
h
m
m hTT
HW,
p.175
#4
§ 4 高斯型 积分 /* Gaussian Quadrature */
ba nk kk xfAdxxfx 0 )()()(?
构造具有 2n+1次代数精度的求积公式将节点 x0 … xn 以及系数 A0 … An 都作为待定系数。
令 f (x) = 1,x,x2,…,x2n+1 代入可求解，得到的公式具有 2n+1 次代数精度。这样的节点称为 Gauss 点,
公式称为 Gauss 型求积公式 。
例,求 的 2 点 Gauss 公式。dxxfx )(1
0?
解,设,应有 3 次代数精度。1
0 1100 )()()( xfAxfAdxxf
x
代入 f (x) = 1,x,x2,x3
3
11
3
009
2
2
11
2
007
2
11005
2
103
2
xAxA
xAxA
xAxA
AA
2776.0
3891.0
2899.0
8212.0
1
0
1
0
A
A
x
x
不是线性方程组，
不易求解。
§ 4 Gaussian Quadrature
证明：,?” x0 … xn 为 Gauss 点,则公式至少有 2n+1 次代数精度。
ba n
k
kk xfAdxxfx
0
)()()(?
对任意次数 不大于 n的多项式 Pm(x),Pm(x) w(x)的次数不大于 2n+1,则代入公式应 精确成立,
n
k
kkmk
b
a m xwxPAdxxwxPx 0 )()()()()(?
0
= 0?
“?” 要证明 x0 … xn 为 Gauss 点，即要证公式对任意次数 不大于 2n+1的多项式 Pm(x) 精确成立，即证明：
n
k
kmk
b
a m xPAdxxPx 0 )()()(?
设 )()()()( xrxqxwxP
m
bababa m dxxrxdxxqxwxdxxPx )()()()()()()(
0
n
k
kk xrA
0
)(
n
k
kmk xPA
0
)(?
x0… xn 为 Gauss 点? 与任意次数不大于 n的多项式 P(x) （ 带权 ） 正交 。
n
k
kxxxw
0
)()(
定理求 Gauss 点? 求 w(x)
§ 4 Gaussian Quadrature
正交多项式族 {?0,?1,…,?n,… }有性质：任意次数不大于 n 的多项式 P(x) 必与?n+1正交。
若取 w(x) 为其中的?n+1，则?n+1的根 就是 Gauss 点。
再解上例,10 1100 )()()( xfAxfAdxxfx
Step 1,构造正交多项式?2
设 cbxxxaxxx 2210 )(,)(,1)(
53a0)(10 dxaxx0),( 10
1
021
1
0
2
10
0))(53(0),(
0)(0),(
dxcbxxxx
dxcbxxx
21
5
9
10
c
b
即：
21
5
9
10)( 2
2 xxx?
§ 4 Gaussian Quadrature
Step 2,求?2 = 0 的 2 个根，即为 Gauss 点 x0,x1
2
21/20)9/10(9/10 2
1;0
x
Step 3,代入 f (x) = 1,x 以求解 A0,A1
解 线性 方程组，
简单。
结果与前一方法相同,2776.0,3891.0,2899.0,8212.0
1010 AAxx
利用此公式计算 的值? 1
0 dxex
x
2555.1 1
0 dxex
x 2899.08212.010 2 7 7 6.03 8 9 1.010 eeeAeA xx
注,构造正交多项式也可以利用 L-S 拟合中介绍过的递推式进行。
§ 4 Gaussian Quadrature
特殊正交多项式族：
① Legendre 多项式族,1)(?x?定义在 [?1,1]上，
k
k
k
kk xdx
d
kxP )1(!2
1)( 2 满足：
lk
lkPP
k
lk
12
2
0),(
xPP 10,1由 有递推 11 )12()1( kkk kPxPkPk
以 Pn+1 的根为节点的求积公式称为 Gauss-Legendre 公式 。
② Chebyshev 多项式族,21 1)( xx定义在 [?1,1]上，
) a rcc o s( co s)( xkxT k
Tn+1 的根为
22 12c o s nkx kk = 0,…,n
以此为节点构造公式
1
1 02 )()(1
1 n
k
kk xfAdxxfx
称为 Gauss-Chebyshev 公式 。
注意到积分端点?1 可能是积分的 奇点,用普通 Newton-Cotes公式在端点会出问题。而 Gauss公式可能避免此问题的发生。
其它公式见教材 p.189
§ 4 Gaussian Quadrature
Gauss 公式的余项：
ba nk kk xfAdxxffR 0 )()(][ /* 设 P为 f的过 x0 … xn的插值多项式 */
ba nk kk xPAdxxf 0 )()(
/*只要 P 的阶数不大于 2n+1，则下一步等式成立 */
dxxPxfdxxPdxxf baba ba )]()([)()(
插值多项式 的余项Q,什么样的 插值多项式 在 x0 … xn 上有 2n+1 阶？
A,Hermite 多项式！ 满足 )()(),()(
kkkk xfxHxfxH
ba dxxHxffR )]()([][
),(,)(
)!22(
)(
)(
)!22(
)(
2
)12(
2
)12(
badxxw
n
f
dxxw
n
f
b
a
n
b
a
x
n
HW,
p.190 #6
