数值分析数值积分

数值分析数值积分第1页/共69页对f()采用不同的近似计算方法，从而得到各种不同的求积公式。第2页/共69页第3页/共69页以上三种方法都是用被积函数值的线性组合来表示积分值。推广，一般地有求积节点求积系数，与被积函数无关像这样，将积分用若干节点上被积函数值的线性组合来表示的数值积分公式称为机械求积公式。求积误差机械型求积公式的构造归结为，确定求积节点xk和求积系数Ak，使在某种意义下精确度较高。总之，要解决三个问题：精确度的度量标准；如何构造具体的求积公式；具体求积公式构造出来后，误差如何估计？第4页/共69页定义：代数精度若某个求积公式对次数m阶的多项式准确成立，而对m+1

阶的多项式不一定准确成立。即对应的误差满足：R[Pk

]=0对任意

k

m阶的多项式成立，且R[Pm+1

]0对某个

m+1阶多项式成立，则称此求积公式的代数精度为

m

。代数精度与误差的关系：代数精度越高，求积误差越小。结论：问题1第5页/共69页第6页/共69页由上面代数精度条件确定求积公式可分两种情形：若事先给定求积节点xk(k=0,…,n),例如被积函数以表的形式给出时xk确定，可令m=n,由上式确定n+1个系数Ak即可----待定系数法和插值法。若xk和Ak都可选择，令m=2n+1，确定xk和法Ak---Gauss法要使求积公式具有m阶代数精度，则它对1,x,…,xm均准确成立，即m+1个方程，2n+2个未知数问题2第7页/共69页Case1---方法1第8页/共69页第9页/共69页§1插值型求积

公式思路利用插值多项式

则积分易算。

在[a,b]上取ax0<x1<…<xn

b，做f的

n

次插值多项式，即得到Ak由决定，与无关。节点

f(x)插值型积分公式/*interpolatoryquadrature*/误差Case1---方法2第10页/共69页§1Newton-CotesFormulae例：对于[a,b]上1次插值，有考察其代数精度。f(x)abf(a)f(b)梯形公式/*trapezoidalrule*/解：逐次检查公式是否精确成立代入P0=1：=代入P1=x：=代入P2=x2：代数精度=1第11页/共69页§1Newton-CotesFormulaeTh1.形如的求积公式至少有n

次代数精度

该公式为插值型（即：）第12页/共69页第13页/共69页当节点等距分布时：令Cotes系数注：Cotes系数仅取决于n

和i，可查表得到。与f(x)及区间[a,b]均无关。§2Newton--Cotes

公式Newton—Cotesformula第14页/共69页§1Newton-CotesFormulaen=1:TrapezoidalRule/*令x=a+th,h=ba,用中值定理*/代数精度=1n=2:Simpson’sRule代数精度=3第15页/共69页n=4:CotesRule,代数精度=5,偶数阶N-C公式具有n+1阶代数精度n=3:Simpson’s3/8-Rule,代数精度=3,第16页/共69页对称节点的系数相同Cotes公式是用不同节点的函数值（高度）的加权平均来近似区间的平均高度注：当n8时，Cotes系数有负，造成公式不稳定，因此常用低阶Cotes公式。第17页/共69页证明：只需证明n为偶数时，N-C公式对f(x)=xn+1的余项R(f)=0即可。因f(n+1)(x)=(n+1)!,由余项公式得Th2.n为偶数时，

N-C公式至少具有n+1阶代数精度。第18页/共69页注：当n

为偶数时，Cotes公式具有n+1阶精度，与n+1阶Cotes公式精度相同，但少计算一个节点上的函数值，因此一般常用偶数阶Cotes公式。第19页/共69页偶数阶N-C公式具有n+1阶代数精度N-C公式具有n阶代数精度余项R=o(h

n+2)第20页/共69页第21页/共69页第22页/共69页第23页/共69页Hint：constructainterpolationpolynomialoforder5,H(x),satisfying

H(a)=f(a),H(b)=f(b),

H(k)((a+b)/2)=f(k)((a+b)/2).第24页/共69页HW:p.151-152#1-6第25页/共69页数值稳定性的一般概念第26页/共69页N-C的稳定性第27页/共69页第28页/共69页§3复合求积

/*CompositeQuadrature*/Haven’twehadenoughformulae?What’supnow?Ohcomeon,youdon’tseriouslyconsiderh=(ba)/2acceptable,doyou?Whycan’tyousimplyrefinethepartitionifyouhavetobesopicky?Don’tyouforgettheoscillatorynatureofhigh-degreepolynomials!Uh-oh高次插值有Runge现象，故采用分段低次插值分段低次合成的Newton-Cotes

复合求积公式。复合梯形公式：在每个上用梯形公式：=

Tn/*中值定理*/第29页/共69页第30页/共69页§2CompositeQuadrature复化Simpson公式：44444=

Sn注：为方便编程，可采用另一记法：令n’=2n为偶数，这时，有第31页/共69页第32页/共69页§2CompositeQuadrature收敛速度与误差估计：定义若一个复化积分公式的误差满足且C0，则称该公式是p

阶收敛的。/*中值定理*/类似的，可得2阶收敛4阶收敛6阶收敛第33页/共69页例1：计算解：其中=3.138988494其中=3.141592502运算量基本相同第34页/共69页Q:给定精度，如何取n?例如：要求，如何判断n=？上例中若要求，则即：取n=409第35页/共69页第36页/共69页第37页/共69页第38页/共69页第39页/共69页第40页/共69页§2CompositeQuadrature第41页/共69页事后误差估计式，可用来判断迭代是否停止。始步长h第42页/共69页第43页/共69页第44页/共69页0.510-2第45页/共69页§4龙贝格积分

/*RombergIntegration*/复化梯形公式算法简单，但精度较差，收敛速度（2阶收敛）较慢，如何提高收敛速度？第46页/共69页第47页/共69页注：按上面规律，可以构造线性组合系数为的新的积分公式，但当m4时，前一个系数接近于1，后一个系数接近于0，这样构造出的新公式与前一个公式结果差别不大，反而增加计算量，因此实际上常做到Romberg公式为止。第48页/共69页第49页/共69页第50页/共69页第51页/共69页第52页/共69页例：计算已知对于=106

须将区间对分9次，得到T512=3.14159202由来计算I

效果是否好些？考察=3.141592502=S4一般有：Romberg序列

Romberg

算法：<?<?<?………………

T1=)0(0T

T8=)3(0TT4=)2(0T

T2=)1(0T

S1=)0(1T

R1=)0(3T

S2=)1(1T

C1=)0(2T

C2=)1(2T

S4=)2(1T第53页/共69页§3RombergIntegration

理查德森外推法/*Richardson’sextrapolation*/利用低阶公式产生高精度的结果。设对于某一h

0，有公式T0(h)

近似计算某一未知值I。由Taylor展开得到：T0(h)I=1h+2h2+3h3+…i与h

无关现将

h对分，得：()()()...)(3232222120+++=-hhhhITaaaQ：如何将公式精度由O(h)提高到O(h2)？...432112)()(23322020---=---hhIhTThaa即：第54页/共69页§5高斯型积分

/*GaussianQuadrature*/构造具有2n+1次代数精度的插值型求积公式令f(x)=1,x,x2,…,x2n+1

代入可求解，得到的公式具有2n+1次代数精度。例：求的2点Gauss公式。解：设，应有3

次代数精度。+101100)()()(xfAxfAdxxfx代入f(x)=1,x,x2,x3不是线性方程组，不易求解。怎么办？先求出Gauss点，则问题转化为前面的问题，或者将方程组变为关于求积系数的线性方程组将节点x0…xn

以及系数A0…An

都作为待定系数。这样的节点称为Gauss点，公式称为Gauss型求积公式。第55页/共69页§4GaussianQuadrature证明：“”x0…xn

为Gauss点,则公式至少有2n+1次代数精度。对任意次数不大于n

的多项式Pm(x)，Pm(x)w(x)的次数不大于2n+1，则代入公式应精确成立：0=0“”

要证明x0…xn为Gauss点，即要证公式对任意次数不大于2n+1

的多项式Pm(x)精确成立，即证明：设0

x0…xn

为Gauss点

与任意次数不大于n

的多项式P(x)（带权）正交。定理求Gauss点

求w(x)第56页/共69页§4GaussianQuadrature正交多项式族{0,1,…,n,…}有性质：任意次数不大于n

的多项式P(x)必与n+1

正交。若取w(x)为其中的n+1，则n+1的根就是Gauss点。再解上例：+101100)()()(xfAxfAdxxfxStep1：构造正交多项式2设cbxxxaxxx++=+==2210)(,)(,1)(jjj53-=a0)(10=+dxaxx0),(10=jj=++-==++=1021102200))(53(0),(0)(0),(dxcbxxxxdxcbxxxjjjj215910=-=cb即：第57页/共69页§4GaussianQuadratureStep2：求2=0

的2个根，即为Gauss点x0，x1Step3：代入f(x)=1,x以求解A0，A1解线性方程组，简单。结果与前一方法相同：

利用此公式计算的值注：构造正交多项式也可以利用L-S

拟合中介绍过的递推式进行（P.102）。第58页/共69页§4GaussianQuadrature特殊正交多项式族：①

Legendre多项式族：1)(xr定义在[1,1]上，满足：由有递推以Pn+1的根为节点的求积公式称为Gauss-Legendre

公式。②

Chebyshev多项式族：211)(xx-=r定义在[1,1