计算方法第4章 - 数值积分

更新时间:2023-10-12 07:48:02 阅读量: 综合文库 文档下载

说明:文章内容仅供预览,部分内容可能不全。下载后的文档,内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的,是否完整无缺。

第四章 数值积分 ( Numerical Integration and Differentiation)

记号: I(f)??f(x)dx——积分, Q(f)??Akf(xk)——求积公式,

ak?0bn误差: E(f)?I(f)?Q(f), ( Quadrature Formula ) 节点:xk, 求积系数:Ak; 4.1 内插求积, Newton-Cotes公式

利用插值多项式p(x):f(x)?p(x)?R(x) 积分 I(f)??p(x)dx??R(x)dx;

aabb例如,插值多项式取Lagrange形式 L(x), 便有

Q(f)??L(x)dx??abbna?lk(x)f(xk)dx??f(xk)?lk(x)dx

k?0k?0anb及

E(f)?I(f)?Q(f)??R(x)dx

ab此类通过在节点xi(i?0,1,?,n),a?x0?x1???xn?b 处满足插值条件的插值多项式导出的求积公式称为内插求积公式。其中:

求积系数 Ak??lk(x)dx,

ab误差 ?R(x)?f(x)?p(x)?f[x0,x1,?,xn,x]?(x)??E(f)??bab1f[x0,x1,?,xn,x]?(x)dx?f(n?1)(?)?(x)dx ?(n?1)!a1f(n?1)(?)?(x)

(n?1)!显然,当函数f(x)为不超过n次的多项式时,必有R(x)?0.由此可见,若函数

f(x)为不超过某次的多项式,便有求积公式的误差E(f)?0,即I(f)?Q(f).

定义:若求积公式对任何不超过m次的多项式pm(x)成立:I(pm)?Q(pm),而对m?1次多项式pm?1(x)等式不再成立I(pm?1)?Q(pm?1),则称该求积公式的代数精度为m (事实上,可对一般的函数与其相应的近似公式定义代数精度).

代数精度——何时没有误差。 本节讨论的公式的节点为等距分布:h?——称为Newton-Cotes 公式

4.1.1 Newton-Cotes 公式

1)n=1,梯形(Trapezoidal)公式 A0?A1?(b?a) ,

2b?a?f(a)?f(b)? T(f)?2

1

1(b?a),xk?a?kh,k?0,1,?,n. n误差:利用积分中值定理(参见附录I,定理I.6:积分中值定理)

E1??f[a,b,x](x?a)(x?b)dx?f[a,b,?]?(x?a)(x?b)dxaabb(b?a)3??f??(?)12

a???b,2)n=2 Simpson公式,

(b?a)a?b[f(a)?4f()?f(b)] S(f)?62(b?a)5(4)误差 E2??f(?)a???b,

2880注:该误差公式不能简单地对仅用此三点的插值多项式的余项公式积分获得。

3)n=4 Cotes 公式

b?a?7f(a)?32f(a?h)?12f(a?2h)?32f(a?3h)?7f(b)? C(f)?908b?a7(6)()f(?)a???b E4??94544.1.2 复化求积公式(Composite Numerical Integration) 梯形公式、Simpson公式的 优点:简单,容易估计误差; 缺点:大区间 误差大 改进:利用积分的可加性

bna?kh?af(x)dx???k?1a?(k?1)hf(x)dx,

1)复化梯形公式

hTn?[(f0?f1)?(f1?f2)???(fn?1?fn)]2 n?1h?[f?a??f(b)?2?f(xk)]2k?1误差:若f??连续,由连续函数的介值定理(见p.329例I.2),及b?a?nh,

h3b?a2ET??[f??(?1)?f??(?2)???f??(?n)]??hf??(?),a???b

1212结论:代数精度 1 ;

计算精度 O(h2)—— 一般,若有误差,误差与步长的关联程度;

E( 若limk?M??,则 E?O(hk) )

h?0h2)复化Simpson公式

x?xkhnSn??[f(xk?1)?4f(k?1)?f(xk)]6k?12x?xkh?[f(a)?f(b)?2?f(xk)?4?f(k?1)];62k?1k?1n?1n

2

误差:若f(4)(x)连续,与复化梯形公式误差的推导一样,有

h5b?a4(4)Es??[f(4)(?1)?f(4)(?2)???f(4)(?n)]??hf(?),a???b;

28802880结论:代数精度 3 计算精度 O(h4);

4.1.3 变步长积分法——步长的选取

通常,对数值积分有一定的误差要求:给定误差界?,使I(f)?Q(f)??; 1) 复化梯形积分法,由误差公式

b?ab?a2I(f)?Tn??()f??(?1)12n

b?ab?a2I(f)?T2n??()f??(?2)122n设f??(x)在[a,b]变化不大,即f??(?1)?f??(?2),将以上两式相除,即得

I(f)?Tn?4

I(f)?T2n由此,得梯形积分的误差估计:

1 I(f)?T2n?(T2n?Tn) ;

31从而,由可计算的值 (T2n?Tn)??3到误差要求。

2) Simpson积分公式,由误差公式

b?ab?a4(4)I(f)?Sn??()f(?1)2880n

b?ab?a4(4)I(f)?S2n??()f(?2)28802n估计

I(f)?T2n??;

由于这是由约等式导出的,因此通常取 T2n?Tn?? 估计数值积分是否达

设f(4)(x)在[a,b]变化不大,即f(4)(?1)?f(4)(?2),将以上两式相除,即得

I(f)?Sn1?16?I(f)?S2n?(S2n?Sn) ;

I(f)?S2n151(S2n?Sn)??15?I(f)?S2n??;

从而,由

3) 变步长 比较步长之比为2:1的复化梯形积分公式:

n?1hTn?[f?a??f(b)?2?f(xk)] ,

2k?1hn?1nx?xk2T2n?[f?a??f(b)?2?f(xk)?2?f(k?1)] ,

22k?1k?1 3

可得(理解):

x?xk1hn1 T2n?Tn??f(k?1)?Tn?(T2n的步长)?(新增点的函数值)

22k?122由此可见,容易通过计算步长对分以后两个不同的复化梯形积分公式计算值的差,估计复化梯形积分公式的误差,由此确定是否达到所要求的精度。

4.1.4 Romberg积分

1由梯形积分的误差估计:I(f)?T2n?(T2n?Tn) ,设想将误差部分“归还”

3给梯形积分,应可以得到积分的一个更精确的近似值。通过计算,有:

I(f)?4T2n?Tn= 3h?n?1nn?1x?x1?h??k2[f?a??f(b)?2f(x)?2f(k?1??)]?[fa?f(b)?2f(x)]?4*???kk?3?222k?1k?1k?1???n?1nx?xkh = [f(a)?f(b)?2?f(xk)?4?f(k?1)]?Sn ;

62k?1k?1注意,梯形积分公式 Tn的精度是O(h2),而通过适当的组合,可得到Simpson公式Sn,它的精度是O(h4),精度得到很大的提高;同时代数精度也得到提高。这种方法并没有增加函数计算量,仅增加了极少的代数运算,就能大大地提高计算精度,这类方法通常被称为“外推法”。

用类似的方法,将Simpson公式适当地组合,可得到精度为O(h6),代数精度为5的 Cotes公式;再将Cotes公式适当地组合,可得到精度为O(h8),代数精度为7的 Romberg公式:

42S2n?SnCn?42?1,43C2n?CnRn? ;

43?1在计算时,可列出如下的表进行计算:

T1T2T4T8T16?S1S2S4S8?C1C2C4?R1R2?

44R2n?Rn一般, Romberg 方法用到 Rn为止,尽管理论上还可以继续进行: 44?1

4

以取得精度更高的公式。但我们可以从两方面衡量该公式的实际意义。

44R2nRn44R2n??R2n,因此可将此式认作是1、将此公式改为:,由于255255255对R2n的修正,而修正值是

Rn;当n较大时,R2n与Rn已是积分的比较准确的255近似值,因此这部分的修正对R2n的进一步改进影响一般比较小;

2、注意到Tn公式的误差:?Th2f??(?),Sn公式的误差:?Sh4f(4)(?),Cn公 式的误差:?Ch6f(6)(?),Rn公式的误差需要用到f的高阶导数,而f的高阶导数的性态是很难确定的(如同插值多项式余项估计),

因此,一般不再继续进行此类外推(Romberg)方法,仅用到 Rn为止。 例 积分:?exdx?1.71828182845904523536028747?

01T1?1.6591409T2?1.7539311S1?1.7188612

T4?1.7272219S2?1.7183188C1?1.7182827T8?1.7205186S4?1.7182824C2?1.7182818R1?1.7182818例: 积分:?T1?0.35777T2?0.43187S1?0.45657T4?0.46030S2?0.46987C1?0.47066 T8?0.47092S4?0.47446C2?0.47477R1?0.47484T16?0.47482S8?0.47612C4?0.47623R2?0.476250.80xdx?0.47703? (效果并不明显)

由此,也可从此体会到,Romberg 方法的效果(准确性),与被积函数的高阶导数的性态有关。

4.1.5 待定系数法

例:求确定以下公式的系数 A,B,C,D 使之具有尽可能高的代数精度,求其代数精度,并求其误差:

?1?1f(x)dx?Af(?1)?Bf(0)?Cf(1)?Df?(0)

解:取Newton插值多项式,先求差商表:

5

证明:1) 在(a,b)中有实根:

反证,设?n?1,若?n(x)?0 (或?n(x)?0)?x?(a,b),则由正交性与

?n(x)?0 (或?n(x)?0)?x?(a,b)性质,有

0?(?n,?0)???(x)?n(x)?0(x)dx???(x)?n(x)dx?0

aabb矛盾,所以至少有一个根。

2)单根:若 x?为重根 ?

?n(x)(x?x?)2为 n?2 次多项式,由正交性

?ba?(x)?n(x)?n(x)(x?x)?2dx?0

?ba?(x)?n(x)?n(x)(x?x?)2dx??ba??(x)??(x)?n??dx?0

?x?x?2矛盾,所以若有根则必为单根。

3)在(a,b)中有且共有n个互异实根: 反证:若仅有l(?n)个根:xi(i?1,2,?,l),即

?n(x)??(x)(x?x1)(x?x2)?(x?xl)

其中 ?(x)?0(或?(x)?0)?x?(a,b)(因为它在此区间中不再有其他根),且是(n?l) 次多项式。将函数

?n(x)(x?x0)(x?x1)?(x?xl)??(x)(x?x0)2(x?x1)2?(x?xl)2

在(a,b)上按权函数?(x)积分,由于(x?x0)(x?x1)?(x?xl)为l(?n)次多项式,因此,按正交性:

b?

a?(x)?n(x)(x?x0)(x?x1)?(x?xl)dx?0

?(x)?(x)(x?x0)2(x?x1)2?(x?xl)2dx?0 (或?0)

另一方面,被积函数在 (a,b) 上保持符号,因此

?ba矛盾。综上,n次正交多项式?n(x)在(a,b) 上恰有n个互异实根。

4.4.2 Gauss型求积公式 1)Gauss 型积分性质

n?1个节点,代数精度达到2n?1的求积公式称为Gauss型求积公式, 这样

的n?1节点称为Gauss点.

问题4:这n?1个Gauss点与正交多项式?n?1(x)的n?1个根关系如何? 定理4.12: 求积公式节点xk:a?x0?x1???xn?b成为Gauss点的充要条件是:xk(k?0,1,?,n) 是(a,b)上以?(x)为权的正交多项式?n?1(x)的n?1个根;

16

证明:? 设xk,k?0,1,?,n为Gauss点,记?n(x)??(x?xk),因代数精

k?0n度:2n?1,故对任何p(x)——不超过n次多项式——

?ba?(x)?n(x)p(x)dx??Aj?n(xj)p(xj)?0

j?0n说明:?n(x)与任何不超过n次多项式正交,??n(x)是第n?1个正交多项式(至多差一个系数?n?1(x)?cn?1?n(x)),即xk,k?0,1,?,n是正交多项式的根。

则任何2n?1?n?1(x)?cn?1?(x?xk),? 设?k(x),k?0,1,?为正交多项式,

k?0n次多项式 q(x)?p(x)?n?1(x)?r(x) —— p(x),r(x):n次多项式;由正交性,

n?1点内插型求积公式的代数精度至少n:

?ba?(x)q(x)dx???(x)p(x)?n?1(x)dx???(x)r(x)dx???(x)r(x)dxaaabbb??Ajr(xj)??Ajq(xj)j?0j?0nn

i.e. 公式代数精度:2n?1——?xk,k?0,1,?,n?为Gauss点。

2) 求积系数及误差

b?n?1(x)?(x)Aj???(x)lj(x)dx???(x)dx???(x)dx

aaa??1(xj)(x?xj)?'(xj)(x?xj)?nbb误差:由于此处考虑一般的Gauss 型积分,因此对代数精度为2n?1的求积公式,求误差公式不是取x2n?2,而是取f(x)???n?1(x)?,此时有Q[??n?1(x)?]?0。

22由于代数精度为2n?1, 一般地有 E(f)?rf(2n?2)(?),因此

?n?1(x)?]?Q[???n?1(x)?]???(x)???n?1(x)?dxI[??a22b2

?n?1(x)?]?r[???n?1(x)?]?E[??22(2n?2)?r?(2n?2)!可知

f(2n?2)(?)b2? E(f)??(x)?(x)dx ??n?1?a(2n?2)!4.4.3 若干典型公式

(各公式的Gauss点,对应的求积系数及误差,请参见教材) (1) Legendre 多项式、积分——权 ?(x)?1,区间 (?1,1),

1dnPn(x)?n[(x2?1)n],n?0,1,2,? n2n!dx或递推公式:

P0(x)?1,P1(x)?x, 2n?1nPn?1(x)?xPn(x)?Pn?1(x),n?1,2,?n?1n?1

17

(2) Laguerre 多项式、积分——权 ?(x)?e?x,区间 (0,??),

nxd?xnLn(x)?e(ex),n?0,1,2,? ndx或递推公式:

L0(x)?1,L1(x)?1?x,Ln?1(x)?(2n?1?x)Ln(x)?xLn?1(x),n?1,2,?(3) Chebychev 多项式、积分——权 ?(x)?2

121?xTn(x)?cos(narccosx),n?0,1,?x?[?1,1]

,区间 (?1,1)

或递推公式:

T0(x)?1,T1(x)?x,Tn?1(x)?2xTn(x)?Tn?1(x),n?1,2,?

算法——利用已有积分公式,作变量代换;

或构造相应的正交多项式,求其根及积分常数。

例:试求以下求积公式的结点xi(i?0,1,?,n)及求积系数Ai(i?0,1,?,n),使公式具有最高代数精度(n?0,1):

?10f(x)dx??Aif(xi)

i?0n解:具有最高代数精度的求积公式必是Gauss型求积公式,对?(x)?1,积分区间(0,1),先求相应的正交多项式:

内积:(f,g)??f(x)g(x)dx

01?0(x)?1

?1(x)?x??0?x? (?0,?0)???1(x)?x?(x?0,?0)

(?0,?0)(x?0,?0)??10dx?1?10xdx?1 21 2?2(x)?(x??1)?1(x)??1?0(x)?[x?111(x?1,?1)??x(x?)2dx?0224 111(?1,?1)??(x?)2dx?02121111??2(x)?(x?)(x?)??x2?x? (可验证其正交性)

2212611n?0,结点: , 插值多项式:p(x)?f()

22

18

(x?1,?1)(?,?)]?1(x)?11?0(x)

(?1,?1)(?0,?0)因此,Gauss 积分公式(即:中矩形公式)为

11 ?f(x)dx?f()

02代数精度为1.

1111 ?,x1??223223x?x0x?x11插值多项式:p(x)? f(x0)?f(x1) 由x1?x0?x0?x1x1?x031113求积系数 A0?3?(x1?x)dx??(x1?x)2?

002211132 A1?3?(x?x0)dx?(x?x0)?

0022n?1,结点:x0?因此,公式当n?1,其Gauss求积公式: 111111f(x)dx?[f(?)?f(?)] ?02223223可验证代数精度?3:

f(x)?1f(x)?xf(x)?x2f(x)?x3f(x)?x41?11111111?[(?)?(?)]?22223223211112112?[(?)?(?)]?3222322311113113?[(?)?(?)]?4222322311114114?[(?)?(?)]522232231 314f(4)(?)1~f(4)(?)122[?(x)]dx?[(x?x)(x?x)]dx 误差:E(f)?101??004!4!f(4)(?)111112?[(x??)(x??)]dx?04!223223f(4)(?)11212?[(x?)?]dx4!?0212f(4)(?)1141121?[(x?)?(x?)?]dx?04!262144

(4)f(?)2152131?[()?()?]4!523?62144f(4)(?)111?[??]4!8072144f(4)(?)?4!?180

19

本例也可用待定系数法解,但n?1时,求解较困难。此外也可利用Legendre积分,经过积分限的变换获得:

对积分?f(x)dx,利用变量代换:x?01因此

t?1 显然,x?[0,1]对应于t?[?1,1], 211t?111f()dt?g(t)dt ?0???1?1222131由Legendre 多项式 P2(t)?t2?,可知结点为:?,通过公式或查

223表可知两求积系数均为1,因此n?1时

1f(x)dx?

?1011t?111111f(x)dx??f()dt??g(t)dt?[g(?)?g()]

2?122?1233t?1t?1)及x? ,可得同样的结果: 22111111)?f(?)] ?f(x)dx?[f(?02223223再利用g(t)?f(

4.5 数值微分

4.5.1 插值公式方法

数值导数与数值积分类似,可由插值公式的导数得到。 先引入一个公式:

df[x0,x1,?,xn,x]?f[x0,x1,?,xn,x,x]: dxf[x0,x1,?,xn,x??x]?f[x0,x1,?,xn,x]df[x0,x1,?,xn,x]?lim证:

?x?0dx?x ?limf[x0,x1,?,xn,x??x,x]?f[x0,x1,?,xn,x,x] ;

?x?0类似的,有二阶导数公式:

d2df[x,x,?,x,x]?f[x0,x1,?,xn,x,x]01n2dxdx

f[x0,x1,?,xn,x??x,x??x]?f[x0,x1,?,xn,x,x]?lim?x?0?x?lim?x?0f[x0,x1,?,xn,x??x,x??x]?f[x0,x1,?,xn,x??x,x]

?xf[x0,x1,?,xn,x??x,x]?f[x0,x1,?,xn,x,x]

?x ?lim?x?0 20

本文来源:https://www.bwwdw.com/article/yplf.html

Top