chap7微分与数值积分

第七 数值微分与数值积7‐1：已测得某地大气压强随高度变化的一组数据高度0压强kgf1000米处的大气压强变化率。注释：已知函数点的值,求在某些点的导数bI[f]b

f(x)dxF(b)Ffxsinx,fxex2x 我们是否可以找简单方法计算下列积分11 11 x ， xe2xdx dx 01

x22 x(2x)2

dx？？x12345y468‐yfx未知只知道部分对应点x12345y468注释：通过函数点的值,求某个区间上的积分积分=已知数据点上函数值的线性§7.1x""x""y""BfABfAC

x各种一阶差商f(x)f(xh)fhf(x)f(x)f(xhf(x)f(xh)f

向前差商公式D向后差商公式D中心差商公式Dff(xh)f(x)hf(x)23f(x)f(x)f(xh)f(x)hf(x)23f(x)f(x)向前:fxfxhfx)hf(

向后:fx

f(x)f(xh)

f(

f(xh)f(xh)

f( 中心:f 二阶差商f''(x)f(xh)2f(x)f (x)f(xh)2f(x)f(xh)2f(4)注意例设f(x=ex分别取步长h102103,104106用向前差商公式和中心差商公式计算f(1.15)hxy向前差商中心差商上表结果说明当步长由102缩小到104时,误差在缩小.进而缩小到106时,由于有效数字丢失,误差反而增大。当步长不太大时,中心差商公式精度优于向前差商公式。以插值多项式的导数作为函数导数的近似,f(xi)n(xiBfABfAC

x插值函数的误差公式f(n1)R(x)f(x)n(x) (n1)!插值型微分的误差公R(x)f(x)(x)

(x)

f(n1)()

(n1)

对任意x[a,] 因(依赖于x)未知,故上式很难估计误差,但若只求某个节点上的导数值,误差可估计.因此,插值型求导公式通常用于求节点处导数的近似值。

R(

)f(

)n(xi)

(n

(xixjji1等距节点常用公式1.两点公1

x0,x1x0f(x)

f(

)(xx1

f(x

(xx0 (

x1

(

x0f(x)L(x)

f(x1)f(x0

i R(x)f(x)L(x)f(0)(xx)hf(

R(x)f(x)L(x)f(1)(xx)hf(

ff(x)f(x1)f(x0)0hR(x)hf(1020f(x)f(x1)f(x0)1hR(x)hf( 21,(a, 2.三点公

xix0

(i二次插值公f(x)f(x

(xx1)(x

f(x

(xx0)(x

f(x

(xx0)(x0

x)(xx

1

x)(xx

2

x)(xx

f(x)(xx1)(xx2)f(x)(xx0)(xx2)f(x)(xx0)(x

插值公式求L(x)f(x)2x

f(x

2x

f(x)2x

因此f(x)L(x)f(x)2x0x1

f(x)2x0x0

f(x)2x0x0

f(x0

h2f

f(x2)

3f(x0)4f(x1)

f(x2f(x)L(x)f(x)2x1x1

f(x)2x1

f(x)2x1x0 hf(x)

f(x)

f(x2)f(x0

f(x)L(x)

f(x)2x2x1

f(x)2x2x0

f(x)2x2x0

f(x)

f(x)

f(x)

f(x0)4f(x1)3f

f(x f(xf(x)

3f(xi)4f(xi1)f(xi2 i f(xi)

f(xi2)4f(xi1)3f(xix i插值公式二次求导L(x)f(x f(x)2f(x

f(x)f(x0)2f(x1)f(x2 i0,1, xix0

(if(xf(x)13f(x)4f(x)f(x0012f(x)1f(x)f(x102f(x)1f(x)4f(x)3f(x2012(x0) (002R(x)2

f(3)( 2R(x)

f(3)( i(a,

(i0,1,ff(x)(x)1f(x)2f(x)f(xi2i012设Sx)fx)的三次样条插值函数,由三次样条f(x)S(且有误差估

xR(k)(x)f(k)(x)S(k)(x)O(h4k (kf(x)f(xh)f(f(x)hf(x)f(xh

f(x)f(x)13f(x)4f(x)f(x0012f(x)1f(x)f(x102f(x)1f(x)4f(x)3f(x2012f''(x)f''(x)f(xh)2f(x)f(x§7.2Newton－Cotesfx)的原函数不存在或计算繁琐,或fx)只有b离散数据点,求I(f) f(x)dx的近似数值fn(x)bI(fn) fn(x)dx作为If)bEn(f) f(x)fn(x) n nbI(f)

f(x)dx

Ln(x)dxlk(x)dxf(xkk 插值型求积公式b

I(f)

f(x)dxAkf(xkkkAkk

lx)dx(k01n)

f无关,称为求积系求积公式的截断误

f(n1)(R(f)R(f)f(x)dxL(x)dx (

(n

xa (k0,1,",n),hbk b

x Ak

lk(x)dx

j

xxjxxa

j t

(1)nk(b

kk

hdt

k!(nk

0(tjkk

(1)nkjnk!(nk)!nn

(tnn

jk(ba)Ckb

n(I(f)

f(x)dx(ba)k

f(xkC(n)

(1)nk k!(nk)!n

(tjjbn阶N-C公式的截断bRn(f)

f(n1)((n

(xxj)dxnn

nh(n1)!nh

f(n1)()(tnnn12Cn12C(nk12123879038790n 梯形公I(f)baI(f)baf(a)4f(ab)f62I(f)f(x)dxbaf(a)fba2梯形公式

f(x)dxhf(x)f(x)b 2b

Simpson公式

f(x)dxhf(x)

4f(x)f(xb 3b

Simpsons3公式8bf(x)dx3hf(x)3f(x)3f(x)f(xb 8

衡量插值型求积公式的精度,用多项式次数作标(一)求定义7.1

f(x)dxAkf(xkbnkbnbn f(x)dxAkf(xkbnk存在某个m1次多项式Pm1x),bb

AkPm1(xknkn则称该求积公式具有m次代数精确度结论

f(x)dxAkf(xkbnkbn

具有m次代数精度n bxldxAxn

(l0,1,"ma k m

m

dx

k

Akxk由n次插值导出的求积公式至少具有n次代数精度bRn(f)b

f(x)dx

nkn

Akf(xk)

b(n n1(bnR(xl)n

(l0,1,",结论2n阶NC公式至少具有2n1次代数精度 设2n1P2n1x

x2n1

"ax

P(2n1)

bb bxanh

nn

(

h2n

-

njn

(t+n-j)

h2n

-

t(t

1)(t222)"(t

n2)对称区间上奇函数积分为 证bf(x)dxb

ba

f(a)4f(ab)fa

ba a 记I(f

f(

I(f)

f(a)4f )f I(x)1(b2

I~

ba(141)ba6ba(a2a2bb)1(b2I I(x2)1(b3a3)

ba

ab

2 I(x)

a

)b

(ba)I(x3)1(b4

baa34(ab)3b3

4a4)4 I(x4

I(x4)1(b5

4ba

44(ab)4b4I(x5

ba(5a44a3b6a2b24ab35b4) 例:确定系数使下列求积公式代数精确度尽可11

f(x)dx

f(1)

f(0)

ffx1xx2时公式等号成立 A1A0A1

A1

A0 A

A f(x)dx

f(1)

f(0)

f

11

fx1xx2精确成立fx)x3时,

A0 A当f(x)x4 左边2右边2

定理7.2fxC(2abbb

f(x)dx

f1b证明：R(f1b

f(x)(xa)(x2bf(b

(xa)(xf()(ba)3 7.3fxC(4ab则SimpsonbR(f) f(x)dxb

baf(a)4f(ab)f(b)

(b 证明

f(4)() f(4)( (a, Simpson公式的代数精度为3寻找fx)插值多项式使得在a,abb点的值等于函数2值且是三次多项式。取fx)的三次Hermite插值多项式H3x)H3(a)f H(ab)f(a H(ab)f(a baf(a)4f(ab)f baH(a)4H(ab)H(b) bH( f(4)( af(x)H3(x) (xa)(x

)2(x bR(f) f(x)dxb

baf(a)4f(ab)f(b)

f(4)(

a f(x)dx

H3(

(xa)(x

)2(x2f(4)(

(xa)(x

a

)2(xb bxabtb f(4)( ba

)1(t1)t(tf(4)((b)5f(4)((b)522253(b

)

f(4)(

f(4 (a f(n1) f(n1)(R(f)ban1(x(n为奇数n(n(n2(n2)! ba(xa2(x(n为偶数其 n1(x)(xx0)(xx1)"(xxn (a, xka(k用 1"，xn

n(k0,1,",

xn1II(f)baf(x)dx(ba)bkf(xk)En(fknb bl(

(k0,1,", ba 次数;但是高次插值有Runge现象,不一定能提复复化求积基本思用分段低次多项式近似被积函数求积分近似复化梯形fx) a

(k0,1,",

yhbyhbnnbf(x)dxb

xk

f(

hf(x)f( k xkk

k0

kf(a)f(b)2

f(xk k bbf(x)dxhf(a)a2f(b)f(x)kkn如果fxC(2)a其截断误差RR(f) (ba)f(2T(a, h

RT(f)

f(x)dx

f(a)f(b)n2 kn2

f(xkh h

f

)]

h2(ba)

f(

(x, k

k

kh2

(ba)f(

(a,fhbn分段Simpson公hbny hx1

h

n

hba

a

(k0,1,",x2x2k

f(x)dx2h[f( 2k

)4f(

2k1)f(x2k44444 44444mbf(x)dxm

f(x)dx2hf

)4f(x )f(x

m hf(a)f(b) m 3 f(x2k f(x2k1 k k baf(x)dxh3mf(a)f(b)kf()kmf(2k)RR(f)sbh (a, h

Rs(f)

f(x)dx

3f(a)f(b)

f(

)4f(x2k1m m

5

kkf(4)(k

k f(4)()

b

k

b

kk

h k(x2k2,x2k

(a,

0

x21问积分区间要等分为多少份才能保证计算结14位有效数字解:0.3e

x0,1

1ex2 0 11042

f(x)ex

f(x)2xef(x)(4x22)ex22f(4)(x)(16x448x212)e

f(x)(8x312x)e

Rs(f

f(4)(

1800

f(4)(

104n421049

n

n8即可1ex2dx0

1 (14e1 38

e1)f(x)(4x22)e

f(x)(8x312x)ex2

fRT(f)2

f(

120 h21104 n n7.3.3逐次分半算法复化梯形公式误差h h (f)

f

[f()

k

k

kbb

f(x)dx

h[22

(b)

注意到区间再次对1h RT(f) [f(b)f(a)]4RT(f2 122 I13I13TnI kks (f)sn

m590km5

f(4)(

h4

k

(4)(kbh4kb180

f(4)(x)dx

(f(b)f(a))注意到区间再次对f(b)f(a)h (f)2

(fSnS I1(I1(SnI 通常采取将区间不断对分的方法即取n取n2m(m0,1,2,"), ba

xa (k0,1,",2m 2m m m

f(a)f(b)2

f(x

(m0,1,2," 2

k 12mf(a(2k1)hm m(m1,2,"k即Tm是前次计算的近似值Tm1的一半与新分点处函数2 2m之和乘以新步长h之和m

Tm称为梯形值序列IR(fIR(f )13)Sm2 Sm2 3f(a)f(b)f(a2kh)mf(a(2khm(m1,2,"2

fxC(4ab2IS R(f,2

)1(

Sm122S S 例:计算椭 4

的周长,使结果具有5位有效数解x2cosy lds4 4sin2cos2d4

3sin2L

则 0

2

|e(I) 813sin213sin2T(12) T1Tf

)

1|TT|

T1T f(3)

8f(8 1|TT| T1T(f()f(3)f(5)f(7)) 13|T4T8| 1

184T842.42211l4T820203sin2 (f(0)4f()

f( ))2. (f(0)4f(

)2f()4f(3)f())

1|SS|

[f(0)4f

)2f )4f( 2f(4

)4f(5)2f(3)4f(7)

f(

1|SS

1104 l4I42.42211梯形公式算法输入ab,f取m1,hba Th[f(a)f 取F0,对k1,22m1,Ff(a(2k1)hT1T 5若|TT|输出I停止m1h TT0转Simpson公式算法输入ab,fx取Ff(a)f(b),Ff(ab), ba(F4F m hb4

3F303

k1,2,"

m1,

f(a(2k1)h)Sh(F2F4F |SS|15ISm1h F2F3 S 转 计算圆周率的一种算法n nsinn 7∵sinxx1x31x51x 7

(1

(1

(1)6 7 逐次分半时粗细近似值的组合 1 2 2 (22

1(1615215

(1)n2

1(42nn3(232

1(1615215

(1)nnn(1)1 2n (2)1 2n 6II(h)a

a

"ahPk 其中01P2" II(h)ah2ah4ah6 问题O(hp1)提高到O(hp2)?

II(h)ah2ah4ah6 II(h/2) 4a216a364

则(2)4 3I[4I(h/2)I(h)]3a24

15a3h6I4I(h/2)I(h)c

ch6 I(h

(h

ch6" 22m (h/2)

2(m1)似Im(h) 22m 则IIm(h) 1对 II(h)a1

ahP2"a 1II(h)aP1 1

aP2

"aPk P1IP1I2k3k2a2

P1)hP2a

P1

"a

P1)hPk(I(h)P1I(h))I (1P1)

a

"a

"(1P1

kI

f( 记 T(k)I(bab

(k0,1,2," 由Euler－MaclaurinIT(k)a(bIT(k)a(ba)2a(ba)4"a(ba)2i012iT(k) 4

则 R(f,T(k))O((ba)4 T(kT(k) m m4mT(k1)T(km4mRf

(k))mm

)2m2m上式Tkm

(m12称为Romberg求积公式m以T (m0,1,2,")作为积分近似值,称为Romberg方法m m停止准则：T(0)T( m 逐行自左向右计分别将求积区间2k等分（逐次对分 T(k

42T(k1)

43T(k1) T01T010 T

T(

42

43 T24T240

T51T5

T632T637 70

818

929

(k)是Simpson序列 2kmmmRomberg方法产生的序列(0当m时收敛于Ifm且收敛速度快于复化Simpson方 一般公 I(f)b

f(x)dxAkf(xkkk求积系Akk

lx)dx(k01n)与f等距节点：N-C公 逐次分半梯形,Simpson 启示：在节点数n固定时适当选取xk}和求积系数Ak可以提高数值积分方法的近似程度 基本思想节点数n固定适当选取xk}和求积系数Ak使求积公式

(x)f(x)dxbnkbn

f(xk具有最高的代数精确度其中x0为权函数,求积系数Ak(k01n)与f无关。 n个节点的求积公式最高代数精确度是多大怎样选取节点与求积系数使求积公式具有最高的代数精确度?bn（一）Gaussbn

a(x)f(x)dx

k

精度,则节点xk}和求积系数Ak}应满足方程A1A2"AnAxA "A A

A

"Axm b b

(x)xldx

(l0,1,…,此2n个未知数的m+1个方程可证,当m12n时有解设 (x)(x

x)2(xx)2"(x

x

b a(x)P2n(x)dxbn AkP2n(xk)k所以n2n定义7.2如果一组节点x1x2,xnab能使求积公 a(x)f(x)dxAkf(xkk具有2n1次代数精度则称这组节点为Gauss点上述求积公式为带权函数(x)的Gauss型求积公式II1f(x)dxf(1)f(133I1f(x)dx5f935)89f(0)59f35)两 公式三 公式例11f(x)dxaf 0.6)bf(0)cf 11中的待定系数ab和c使其代数精度尽量高,并,1解:If

f(x)dx,I(f)af 0.6)bf(0)cf 取fx1x,x2使I(fI(f1I(x)

1dx I(1)ab111xdx I(x) 0.6a11

I(x2)

x2dx2 I(x2)0.6a0.6c要使公式具有尽可能高的代数精度ababc

ac9 a

c 2

b0.6a0.6c ac5b8时求积公式的代数精度最高 1f(x)dx5f90.6)8f(0)5f990.6I(x3)

x3dx

5

0.6)3

I(x I(x4)

x4dx5

~(x4)

59

0.62)0.40.60.6 I(x)1x

I(x) 0.6)5 0.6

(6II(x6)1x6dx ~x)5(0.630.63)2.16(6I 4例: 求积公式计算4

xe2xdx解：通过坐标变换化为[-1，11令x2t14I4

xe2xdx

(4t4)e4t

f(t(4t4)e4t4两 公式

4

4I1f(t)dtf

)f

)(4

3(4 3I1f(x)I1f(x)dxf(1)f(133

三 公式

4xe2xdx0I

f(t)dt5f 0.6)8f(0)5f 115(49

0.6)e4

8(4)e45(4

II1f(x)dx5f953)8f(0)5f9935)5 )8 )5 注释:Gauss型求积公式由代数精确度描述 最高阶,达到2n-求积节点和系数有特殊的选取 或者（二）Gauss 确定x1,x2,A1,A2,使得求积公1

f(x)dxA1f(x1)

f(x2具有最高代数精确度。解:Gauss求积公式可以达到3次代数精确度.取fx)为任意三次多项式,利用多项式除法,fx)可表成如下形式f(x)(a0a1x)(xx1)(xx2)(b0

f(x)dx

a1x)(xx1)(xx2)dx

A1f(x1)A2f(x21(011)如果求积公式具有3次代数精确度，则11(b0b1x)dxA1(b0b1x1)A2(b0b1x21 1(a0a1x)(xx1)(xx2)dx对任意实数a0,a1都成立1 (xx)(xx)dx1 x(xx1)(xx2)dx即 22xx 22

(x1x2)3由此解出x= 13 再取两个特殊的f,求得A1和A2:比如取f11f21 dx2A1

331- 33

xdx0 其解为A1A21。故所求求积公式为1f(x)dxf1f13333

例如f(x)xx2xx 则f(x)dx

1x-x2x-

2dx11而

f(x1)f(x2)b求积公式b

(x)f(x)dx

nkn

fxk)nn(x)(xxknk与任意至多n1次多项式qx)在[ab]上关于权函数x)正交即b(x)q(

(x)dx 若xk}为Gaussba(x)f(x)dxb

nkn

Akf(xk为带权函数x)的Gauss型求积公式n则对任意至多2n1Px)nba(x)P(x)dxb

k

AkP(xk若qx)n1次多项式bn则qx)n(x)至多2n1次多项式bn

(x)q(x)n(x)dx

k

Akq(

)n(

)

(x)(xxkkf(xxl(l01n1)A1A2"AnAxAx"A

(x)xldxbab

lnAnA

A

A

"

Axn1

b

b(x)xldx

n

(k1,2,",ka(x)f(x)dxk

k

Akf(

代数精度至少为n设fx)为任意次数不超过2n1则存在次数不超过的n1多项式qx),r(x)使f(x)q(x)n(x)r(bb(x)f(x)dxb

(x)q(

(x)dx

b(x)r(a Akr(xk)Akf(xk bkb

k即求积公式

(x)f(x)dxk

Akfxk)对次数不超2n1次的多项式精确成 求积公式

(x)f(x)dx

k

f(

)为Gauss公式,{xk}是Gauss点 证回定理6.3设nx)(n0,12是最高次项系数不为零的n条件：对任意n1次多项式qx都有bb

(x)q(

(x)dx (n1,2,"求求Gauss knk[ab]上关于权函数x)的正交项式系中的n次多项式的n个实根A b(x)l( A b(

n(

(xx)(x因为lkx)

nnj

xxxkx

是n1次多项式,所jk b

n( a(x)lk(x)dx

(

(xx)(x)dx

Allk(xl)注:

bk(x)l2(bk

lk l2(k

xxj

jj

xkxjnn b(x)l2(x)dx a

l

Al2(x) 结论f(x)C(2n)a,b,则(a, 使bR(f)b

(x)f(x)dx

nkn

Akf(xkf(2n)((2n)!

b(x)2(x 证明：构造过节点x1,x2xn上f(x)的2n-1阶f2n( (x) (a, R(f)

(x)f(x)dxk

Akf(xk

a(x)f(x)dx

k

AkH(xk b(x)f(x)dxb(x)H(

f(2n)(x)bb

(x)2(nf(2n)()bn

n(n

2(

(a, fxC(ab),Gaussn时bb

(x)l(x)dxbkabk

kbakb

nkknk

(nnkn

f(xk)

nkn

Ak

(x

nkkknkk

f(x)f(x k

Akk1kkk

f(xk)

(xbb

(x)dx1k

f(xk)

(x

[1,1上带权函数x1Gauss11

f(x)dxAkf(xknnGauss点xk}n次Legendre多项nPx) x21)nn个零点n 2nn!dxn 求积系

1,2,",n) (1x2)P

(2n

b求积分b

f(

做变

xbatba 11b f(x)dxb

b2

f(batba kbak

Af(ba

ba k

II1f(x)dxf(1)f(133I f(x)dx5f135)89f(0)59f935)11的Gauss型求积公1f(x1xndxAkf(xkkGauss点 cos2k1(k1,2,",n)为n次Chebyshev多项kTnxcos(narccosx)n个零

(k1,2,", f(

2k1x1x

dx

nkn

f(cos 插值余:R(f) f插值余

[0)上带权函数x)ex的Gauss型求积公 f(x)dxAkf(