数值分析题

按定义，上述各数具有5位有效数字的近似数分别是 注意

187.93, 0.037856, 8.0000, 2.7183

x8.000033的5位有效数字的近似数是8.0000而不是8，因

2

重力常数

为8只有1位有效数字。

g，如果以m/s2为单位，

g0.980101m/s2；若以km/s2为单位，

22g0.98010km/s，它们都具有3位有效数字，因为按第一

例种写法

112g9.80101013，

22按第二种写法

115g0.00980101023，

22他们虽然写法不同，但都具有3位有效数字。至于绝对误差限，由于单位不同结果也不同，*11122*10m/s,2105km/s2，而相对误差都是22r*0.005/9.800.000005/0.00980。

例3 要使设取知a20的近似值的相对误差限小于0.1%，要取几位有效数字？

110n1。由于204.42a1，

n位有效数字，由定理1.1，r*14，故只要取n4，就有

r*0.1251031030.1%

即只要对表得20的近似值取4位有效数字，其相对误差就小于0.1%。此时由开方204.472。

例4 设

yxn，求y的相对误差与

x的相对误差之间的关系。

x的相对误差是x的

解 由式(1-9)得

er(y)d(lnxn)nd(lnx)ner(x)

所以

xn的相对误差是x的相对误差的n倍，特别地，相对误差的一半。

例5 设

x0，x的相对误差为，求lnx的绝对误差。

e(x)，即e(x)x，所以

er(x)xe(x)e(lnx)d(lnx)er(x)

x解 由于l100m,宽d的值为

*d*80m,已知|ll|0.2m，|dd*|0.1m，试求面积sld例6 已测得某场地长为

的值为

的绝对误差与相对误差。

解 因sl*ssld,d,l,由(1-5)知

ld*其中(s)*ls**d80m,()l110m，而

d*ss**s||l||d*,

ld*l*0.2m,d*0.1m,

于是绝对误差限

*s800.21100.127m相对误差限

*2

27rs***0.31%

|s|ld8800

例7 计算ns*(s*)Ie1由分部积分可得计算n的递推公式

I10xnexdx(n0,1,)并估计误差。

若计算出0，代入(1-11)，可逐次求出1算

IIn1nIn1 (n1,2,) (1-11)

11x1Ieedx1e.00I,I2,的值。要算出0就要先计

Ie1，若用泰勒多项式展开部分和

(1)2e1(1)2!1(1)k， k!并取

k7，用

14位小数计算，则得

11R7e0.3679104。计算过程中小数点后第5位的数

8!4字按四舍五入原则舍入，由此产生的舍入误差这里先不讨论。当初始值取为

e10.3679，截断误差

I00.6321I0时，用(1-11)递推的计算公式为

I0（A）0.6321; (n1,2,)。 In1nIn1计算结果见表1-1的n列。用0近似0产生的误差

始误差，它对后面计算结果是有影响的。

表1-1

IIIE0I0I0就是初

I*n 0 1 2 3 4 In (用(A)算) 0.6321 0.3679 0.2642 0.2074 0.1704 I* n 5 6 7 8 9 In (用(A)算) 0.1480 0.1120 0.2160 -0.7280 7.552 (用(B)算) 0.6321 0.3679 0.2643 0.2073 0.1708 (用(B)算) 0.1455 0.1268 0.1121 0.1035 0.0684 从表中可以看到8出现负值，这与一切n得

II0相矛盾。实际上，有积分估值

(1-12)

1e11xe(mine)xndxIn00x1n1e(maxe)xndx0x101x11n1n较大时，用In近似In显然是不正确的。这里的计算公式与每步计

算都是正确的，那么是什么原因使计算结果错误呢？主要就是初值I0有误差

因此，当

E0I0I0，由此引起以后各步计算的误差EnInIn满足关系

EnnEn1 (n1,2,)

容易推得

En(1)nn!E0，

由此看出：误差

E0导致第n步的误差扩大n!倍，当n较大时，误差将淹

没真值，因此用n近似In显然是不正确的，这种递推公式不宜采用。例如，

In8，若E0I1104，则 2E88!E02

这就说明8完全不能近似8了。它表明公式(A)是数值不稳定的。

现在换一种计算方法。由(1-12)取

In9，有

111e我们粗略取I(，然后将公式(1-11)倒过来算，)0.0684I9921010e11

I91010**即由9算出8II,I,*7,I*0公式为

*I90.0684（B） n(9,8,*1*I(1In)n1n,1);

4**列。我们发现I0与I0的误差不超过10。记In1******EnInIn，则E0En，E0比En缩小了n!倍，因此，

n!**尽管E9较大，但由于误差逐步缩小，故可用In近似In。反之，当用方案(A)计算时，尽管初值I0相当准确，由于误差传播是逐步扩大的，因而计算结果不计算结果见表1-1的

可靠。此例说明，数值不稳定的算法是不能使用的。

例8 线性方程组

0.00001x1x21,  2x1x22.的准确解为

200000199998x10.50000125,x20.999995

399999199999现在四位浮点十进制数（仿机器实际计算）下用消去法求解，上述方程写成

411100.1000x100.1000x100.1000, 11111100.2000x100.1000x100.2000.1214若用（100.1000)除第一方程减第二方程，则出现用小的数除

2大的数，得到

411100.1000x100.1000x100.1000, 1166100.2000x100.2000.2由此解出

x10, x21010.10001,

显然严重失真。

若反过来用第二个方程消去第一个方程中含1的项，则避免了大数被小数除，得到

66100.1000x100.1000,2 111100.2000x100.1000x100.2000.12x由此求的相当好的近似解

x10.5000,x21010.1000。 2例9 求x16x10的根。

63， 解 x18* x28630.8100.794100.00610x2x*2只有一位有效数字。若改用

1x28630.627101

86315.941具有3位有效数字。

例10 计算A107(1cos2)（用四位数学用表）。

0.9994，直接计算

773A10(1cos2)10(10.9994)610。

x只有一位有效数字，若利用1cosx2sin2，则 2由于cos2A10(1cos2)2(sin1)106.1310具有三位有效数字（这里sin10.0175）。

和

7273

此例说明，可通过改变计算公式避免或减少有效数字的损失。类似的如果1xx2很接近时，则

x1

lgx1lgx2lgx2用右边算式有效数字就不损失。当

x很大时，

x1x都用右边算式代替左边。一般情况，当fxfx*时，就用泰勒展开

fx****2

fxfxfxxxxx2取右端的有限项近似左端。如果无法改变算式，则采取增加有效位数进行运算；在计算机上则采用双倍字长进行运算，但是增加计算时间和多占内存单位。

例11 在五位十进制计算机上，计算

x1x1

A52492i，其中0.1i0.9

i1把运算的数写成规格化形式

1000A0.5249210i

5由于计算机内计算时要对阶，若取

1000i1i0.9，对阶时

i0.000009105，在五位的计算机中表示为机器零，因此

A5.2492100.000009100.000009105555

5.249210（符号

表示机器中相等）。结果显然是不可靠的。这是由于运算中大数“吃掉”

小数

i0.9造成的。

299例12 求二次方程x(101)x100的根。

9解 利用因式分解容易求出此方程的两个根为x110,x21，但若

101(101)410x29929

用求根公式，则得

若用8位小数的计算机运算，由于对阶，有

1091109(101)410109这样求得x110,x20，结果显然是错误的。为避免这种情形出现，

也可采用改变计算公式的方法。如将式

9299

1091(1091)24109x22改变成

x2则有

2109

1091(1091)24109210x291

910109此结果是精确的。

第二章 插值法

sin0.320.314567,sin0.34

0.333487,sin0.360.352274用线性插值及拋物线插值计算sin0.3367的值，并估计截断误差。

例

1

已

给

解 由题意取

x00.32,y00.314567,x10.34,

y10.333487,x20.36,y20.352274

用线性插值计算，取x0.32及x0.34，由公式(2-4)得

01sin0.3367L1(0.3367)0.3367x00.3367x1 y0y1x0x1x1x00.330365其截断误差由(2-9)得

其中

M2R1(x)(xx0)(xx1)

2M2maxf(x)。因f(x)sinxx0xx1x0xx1，可取

M2maxsinxsinx10.3335，有

R1(0.3367)sin0.3367L1(0.3367)10.33350.01670.0033

250.9210用抛物插值计算sin0.3367。由公式(2-4)得

(xx1)(xx2) L2(x)y0(x0x1)(x0x2)(xx0)(xx2)(xx0)(xx1) y1y2(x1x0)(x1x2)(x2x0)(x2x1)有

sin(0.3367)L2(0.3367)0.330374

这个结果与6位有效数字的正弦函数表完全一样，这说明用二次插值精度已相当高了.其截断误差限由(2-9)得

M3R2(x)(xx0)(xx1)(xx2),

6f(x)cosx00.828，于是 其中Mmax3x0xx2R2(0.3367)sin0.3367L2(0.3367)0.17810例2.1 设f(x)x，并已知

6

x f(x) 解 构造均差表如下

2.0 1.414214 2.1 1.449138 2.2 1.483240 试用二次Newton插值多项式

N2(x)计算f(2.15)的近似值，并讨论其误差。

一阶均差

0.34924 0.34102 二阶均差

xk2.0 2.1 2.2 f(xk)1.414214 1.449138 1.483240 0.04110 利用Newton插值公式(2-12)有

N2(x)1.4142140.34924(x2.0)0.04110(x2.0)(x2.1)取x2.15，得N(2.15)1.466292。

2由于f在区间[2.0,2.2]上充分光滑，因此可以利用误差估计公式(2-11)，注意

3(3)(3)到f(x),maxf(x)0.06629。

2(8xx)2.0x2.2从而得到

0.001maxf(x)N2(x)0.066292.0x2.2 430.552417103f(2.15)的真值为1.466288，因此得出

5R(2.15)0.410。由此看出，在较小区间上用式(1.10)，可得到一个

(3)较好估计。

f(x)的函数表2-2，求4次牛顿插值多项式，并由此近似计算

f(0.596)。

例2 给出

首先根据给定函数表造出均差表。

表 2-2 0.40 0.55 0.65 0.80 0.90 1.05 三阶 0.19733 0.21300 0.22863 四阶 0.03134 0.03126 五阶 -0.00012 xkf(xk)0.41075 0.57815 0.69675 0.88811 1.02652 1.25382 一阶 1.11600 1.18600 1.27573 1.38410 1.51533 二阶 0.28000 0.35893 0.43348 0.52483 N4(x)作近似即可。

N4(x)0.410751.116(x0.4)从均差表看到4阶均差近似常数，故取4次插值多项式

0.28(x0.4)(x0.55)0.19733(x0.4)(x0.55)(x0.65) 0.03134(x0.4)(x0.55)(x0.65)(x0.8)于是

f(0.596)N4(0.596)0.63192,

截断误差

R4(x)f[x0,这说明截断误差很下，可忽略不计。

,x5]5(0.596)9

3.6310此例的截断误差估计中，5

f[x,x0,,x4]用

f[x0,,x5]0.00012近似。另一种方法是取x0.596，

,x4]的近似值，从由f(0.596)0.63192，可求得f[x,x,0而可得R(x)的近似。

4例2.4 已知ysinx的函数表如下，分别用牛顿向前、向后插值公式求

阶均差

sin0.57891的近似值。

x 0.4 sinx 0.38942 解 取

0.5 0.47943 0.6 0.56464 0.7 0.64422 h0.1。按表2.4计算得

yi 0.38942 0.47943 0.56464 0.64422 1 0.09001 0.08521 0.07958 x00.4,x10.5,x20.6.x30.7，有

一阶差分 二阶差分 -0.00480 -0.00563 三阶差分 -0.00083 1 t 1t(t1) 2t(t1)(t2) 3! t t(t1) 2t(t1)(t2) 3!Newton向前插值公式为

N3(x0th)0.389420.09001t

xx00.578910.4将t1.7891代入上式得

h0.1110.00480t(t1)0.00083t(t1)(t2) 26sin0.57891N3(0.57891)0.389420.090011.7891110.004801.78910.7891260.000831.78910.78910.2109

0.54711由式(2-26)，误差为

(0.1)R(0.57891)1.78910.78914!(0.2109)(1.2109)sin2106Newton向后插值公式为

4

N3(x3th)0.644220.07958t

0.005630.00083t(t1)t(t1)(t2)23!xx3将t1.2109代入上式得

h1sin0.578910.644220.079581.21092!0.00563(1.2109)(0.2109)查表可得

sin0.578910.5471118。

10.00083(1.2109)3!(0.2109)0.78910.54711

如果取

x00.4,x10.5,x20.6，用二阶Newton向后插值公式，

则得

1N3(x2th)0.564640.08521t0.00480t(t1)2

xx20.2109代入上式，得 将thsin0.578910.564640.08521(0.2109)10.00480(0.2109)0.7891

20.54707其误差为

(0.1)3R(0.57891)0.21090.78911.78913!cos0.510

例

2.5

给

出

4f(x)cosx在xkkh，

k0,1,,6,h0.1处的函数值，试用4次等距节点插值公式计算f(0.048)及f(0.566)的近似值并估计误差。

解 构造差分表2.5。用牛顿向前插值公式(2-25)计算f(0.048)的近似值，取

x00.48，用表2.5上半部差分，x0.048，h0.1，th得

表 2.5

f(xk) f(f) 2f(2f) 3f(3f) 4f(4f) 5f(5f) 1.00000 0.99500 0.98007 0.95534 0.92106 0.87758 0.82534 -0.00500 -0.01493 -0.02473 -0.03428 -0.04348 -0.05224 -0.00993 -0.00980 -0.00955 -0.00920 -0.00876 0.00013 0.00025 0.00035 0.00044 0.00012 0.00010 0.00009 -0.00002 -0.00001 N4(0.048)1.000000.48(0.00500)(0.48)(0.481)1(0.00993)23!(0.48)(0.481)(0.482)(0.00013)1(0.48)(0.481)(0.482)4!(0.483)(0.00012)0.99885cos0.048误差估计由(2-26)可得

M5R4(0.048)t(t1)(t2)(t3)(t4)h5 5!1.5845107,其中Msin0.60.565.

5f(0.566)。用牛顿向后插值公式(4.12) 计算

xx6x0.566,x60.6,t0.34，用差分表2.5中下

h半部差分，得

N4(0.566)0.825340.340.052240.660.008760.000090.000441.662.6626240.84405于是

cos0.5660.84405，误差估计由(2-28)得

M5R4(0.566)t(t1)(t2)(t3)(t4)h5 5!1.7064107,M50.565.

2.6

求

满

足

例

其中

P(xj)f(xj)(j0,1,2)及

P(x1)f(x1)的插值多项式及余项表达式。

解 按插值条件，所求P(x)是一个次数不高于3的多项式，它的曲线过点(x0,f(x0))，(x1,f(x1)),(x2,f(x2))，故可设

P(x)N2(x)C(xx0)(xx1)(xx2)

其中

N2(x)f[x0]f[x0,x1](xx0)f[x0,x1,x2](xx0)(xx1)

C为待定常数。

由条件P(x)f(x)可确定常数C，通过计算可得

11f(x1)f[x0,x1](xx0)(x1x0)f[x0,x1,x2]C(x1x0)(x1x2)P(x)与f(x)的误差函数为

R(x)f(x)P(x)

R(xi)f(xi)P(xi)0 (i0,1,2)以及

R(x1)f(x1)p(x1)0，故可设

2R(x)k(x)(xx0)(xx1)(xx2)

其中k(x)为待定函数。为求k(x)，引进辅助函数

2(t)f(t)P(t)k(x)(tx0)(tx1)(tx2) 显然(x)0 (j0,1,2)。且(x)0，(x)0，故

i1(4)(t)在(a.b)内有五个零点（二重根算两个）。反复应用罗尔定理，得(t)在(a.b)内至少有一个零点，故

(4)()f(4)()4!k(x)0

由

于

由此可得余项表达式为

()2R(x)(xx0)(xx1)(xx2)

4!式中位于x,x,x和x所界定的范围内。

012另一个典型例子是两点三次Hermite插值，插值节点为x及xkk1，插值多项式为H3(x)，满足条件

fH3(xk)yk,H3(xk1)yk1, (4.6)

(xk)mk,H3(xk1)mk1H3采用构造基函数方法。假设我们能够构造出基函数

(4)j(x),j(x),jk,k1，满足条件

j(xi)ij,j(xi)0那么多项式

,i,jk,k1 (4.7) j(xi)ijj(xi)0,H3(x)k(x)ykk1(x)yk1k(x)mkk1(x)mk1

便是满足条件(4.6)的插值多项式。

余下的问题就是如何构造出插值基函数

j,j,jk,k1。由于

k(x)在xk1处函数值与导数值均为0，故它应含因子(xxk1)以设

2，因此可

xxk1k(x)[a(xxk)b]xkxk1由条件(4.7)知(x)1,(x)0，故

kkkkb1 2a0xkxk1由此有

2

xxkxxk1k(x)12xk1xkxkxk1同理

2

xxk1xxkk1(x)12xk1xkxk1xk2xxk1k(x)xxkxkxk12

xxkk1(x)xxk1xk1xk因此节点三次Hermite插值多项式为

2

xxkxxk1H3(x)12ykxk1xkxkxk1xxk1xxk12yk1xkxk1xk1xkxxk1(xxk)mkxkxk1xxk(xxk1)mk1xk1xk其插值余项为

2222

R(x)f(x)H3(x)例2.7 已知函数表

f

(4)()22(xx0)(xx1)4!j xj yj

求满足边界条件数。

0 27.7 4.1

1 28 4.3

2 29 4.1

3 30 3.0

y(27.7)3.0,y(30)4.0的三次样条插值函

解 我们采用三弯矩方程求解。由边界条件，所适用的方程组为

21M0d02Md1111 (2-58)

222M2d212M3d3y3y26y1y06),d3(y3其中d(y0)。

0h0h0h2h2计算果见下表

hj(j0,1,2),j,j (j1,2),dj(j0,1,2)的结

j 0 1 hj 0.30 1 j j dj 46.666 2 3 1 3 131 2 10 4.0002 1312.70000 217.4 将以上数据代入方程组(2-58)，解得

M023.531,M10.395,M20.830,M39.115

将此解代入以弯矩为参数的样条分段表达式，得

13.07278(x28)314.843229(x28)30.21944(x27.7)14.31358(x27.7)x[27.7,28]30.06583(29x)4.23417(29x)3S(x)0.13833(x28)3.96167(x28)x[28,29]0.13833(30x)33.96167(30x)31.51917(x29)4.51917(x29)x[29,30] 例2.8 求满足下面函数表所给出的插值条件的自然样条函数，并算出

f(3),f(4.5)的近似值：

j xj yj 0 1 1 1 2 3 2 4 4 3 5 2 解 我们采用三弯矩方程求解。由边界条件，所适用的方程组为

其中

d02f00,d32f30。造差商表如下：

200d02Md1111 (2-59)

222M2d2020d3jx f(x) f[x,x] f[x,x,x] jjj1jj1jj10 1 2 1 2 4 1 3 4 2 12 12 3 5 2

2

56

h0x1x01,h1x2x12,h2x3x21h0h11122

1,2h0h1123h1h22131d16f[x0,x1,x2]6()3,

25d26f[x1,x2,x3]6()5

6代入方程组(2-59)，解得

39M1,M2

44所以

117318(x1)8(x1), x[1,2)3172S(x)3(x2)(x2)(x2)3,x[2,4)

884935234(x4)(x4)(x4),x[4,5]488经计算得

17f(3)S(3)4.25,f(4.5)S(4.5)3.1406

4例2.9 已知函数表如下：

xi f(xi) 求

0.1 0.70010 0.2 0.40160 0.3 0.10810 0.4 -0.17440 0.5 -0.43750 f(x)在0.3与0.4之间的根x。

解 按距

f(x)0由近及远的顺序排列节点yk，并作如下均差表 yk 0.10810 k 0 1 2 3 4 xk 0.3 0.4 0.2 0.5 0.1 -0.35398 -0.34722 -0.35753 -0.35162 0.023032 0.039163 0.019794 -0.029564 -0.022148 0.012527 -0.17440 0.40160 -0.43750 0.70010 N2(y)x0f[y0,y1](yy0)f[y0,y1,y2](yy0)(yy1)于是

N2(y)0.30.35398(0.10810)0.023032(0.10810)0.17440

0.33831函数在0.3与0.4之间的零点为x0.33831。

注：1．当采用低阶插值函数时，按距f(x)0由近及远的顺序排列节点yk，可使插值误差极小化。

2．若选择所有的节点进行插值，由插值多项式的唯一性，与节点的排列顺序无关。

第三章 函数逼近与曲线拟合

例1

R与C的内积.设x,yR,x(x1,,xn)y(y1,,yn)T，则内积可定义为

ni1nnnT，

(x,y)xiyi (1.4)

由此导出向量2-范数为

x若给定实数

2(x,x)2xi i1ni0 (i1,,n)，称{i}为权系数，则在Rn上可

定义加权内积为

(x,y)ixiyi

i1n (1.5)

相应的范数为

x22xii i1n不难验证(1.5)给出的

(x,y)

满足内积定义3.2的条件.当

i1 (i1,,n)时，(1.5)就是(1.4).

n如果x,yC，带权内积定义为

(x,y)ixiyi

其中

nC[a,b]上定义带权的内积，为此，我们先给出权函数的定义.

例2 C[a,b]上的内积.设f(x),g(x)C[a,b]，(x)是[a,b]上给定的权函数，则可定义内积

也可以在

i仍为正实数序列，yi为yi的共轭.

i1(f(x),g(x))(x)f(x)g(x)dx

ab容易验证它满足内积定义的四条性质，由此内积导出的范数为

bf(x)2(f(x),f(x))(x)f2(x)dx a分别称为带权(x)的内积和范数，特别常用的是(x)1的情形，即

1212(f(x),g(x))f(x)g(x)dx

abb2f(x)2f(x)dxa12

2 曲线拟合的最小二乘法

例1 给定一组数据如下：

i xi yi 求

1 2 1.1 2 4 2.8 3 6 4.9 4 8 7.2 x,y的函数关系。

解 先作草图。如图2所示，这些点的分布接近一条直线，因此可设想

y为x的一次函数。设

ya1xa0 (2.1)

从图2不难看出，无论0选取0

a,a1取何值，直线都不可能同时过全部数据点。怎样

a,a1，才能使直线(2.1)“最好”地反映数据点的基本趋势？首先要建立

好的标准。

8 6 4 2 4 6 图2

a,a1已经确定，yi*a1xia0 (i1,2,3,4)为由近似函数求得的近似值，它与观测值yi之差

*iyiyiyia1xia0 (i1,2,3,4)

假设0称为残差。显然，残差的大小可作为衡量近似函数好坏的标准。常用的准则有以下三种：

(1)使残差的绝对值之和最小，即(2)使残差的最大绝对值最小，即

miniii； ；

minmaxi(3)使残差的平方和最小，即

mini2。

i准则(1)的提出很自然，也合理，但实际适用不方便。按准则(2)来求近似函数的方法称为函数的最佳一致逼近。按准则(3)确定参数，求得近似函数的方法称为最佳平方逼近，也称曲线拟合（或数据拟合）的最小二乘法。它的计算比较简便，是实践中常用的一种函数逼近方法。

例2 求数据表 i xi yi 1 -1 -0.2209 6 0.25 2.5645

2 -0.75 0.3295 7 0.5 3.1334

3 -0.5 0.8826 8 0.75 3.7601

24 -0.25 1.4329 9 1 4.2836

5 0 2.0003

的最小二乘二次拟合多项式。 解 二次拟合多项式为2组(2.3)，可得

P(x)a0a1xa2x，将数据代入正则方程

9a003.75a218.172403.75a108.4842 3.75a02.7656a7.617302其解为0a2.0034,a12.2625,a20.0378，所以此数据组

的最小二乘二次拟合多项式为

P2(x)2.00342.2625x0.0378x例3 设一发射源的发射强度公式形如如下表

2

II0e0.5 1.34 t，现测得

I与t的数据

0.8 0.56 ti Ii 0.2 3.16 0.3 2.38 0.4 1.75 0.6 1.00 0.7 0.74 解 先求数据表

ti lnIi 0.2 1.1506 0.3 0.8671 0.4 0.5596 0.5 0.2927 0.6 0.0000 0.7 -0.3011 0.8 -0.5798 的最小拟合直线。将此表数据代入正则方程组(2.3)，可得

7a03.5a11.9891 3.5a02.03a10.1858其解为a01.73,a12.89。所以

I0e5.64,a12.89

2.89t发射强度公式近似为I5.64e。

例4 在某化学反应里，根据实验所得生成物的浓度与时间关系如下表，要求浓度

a0f(t)与时间t的拟合曲线yy(t)。

1 2 3 4 时间t（分） 34.0 5.4 8.0 3.8 浓度y10 10

10.2

11 10.32

12 10.42

13 10.5

5 9.22 14 10.55

6 9.5 15 10.58

7 9.7 8 9.8 16 10.60

9 10.0

解 将数据描在坐标纸上，如图3。我们看到开始时浓度增加较快，后来增长逐渐减弱，当

于某个常数，故有一水平渐近线。另外，当合曲线

t时y趣

10 8 6 4 2 y103 t0时，反应还未开始，浓

yy(t)是双曲型

ty

a1a0t度应为零。根据这些特点，可设想拟

2 4 6 8 10 12 14 16 t

它与给定数据的规律大致符合。上述模型是非线性参数问题，可以通过变量的变换

图3

11z,x

yt变为线性参数的数学模型

(xi,zi),i1,2,za0a1x拟合数据

,16。其中xi,zi分别由原始数据ti,yi根据变换

公式计算出来。我们建立相应的法方程组

316a3.38073a1.83721001 33.38073a01.58435a10.5288610解此方程组得

a080.6621,a1161.6822

从而得拟合曲线

tyy1(t)

80.6621t161.6822求数据组的最小二乘拟合函数的步骤： （1）由给定数据确定近似函数的表达式，一般可通过描点观察或经验估计得到； （2）按最小二乘原则确定表达式的参数，即由正则方程组，求解得参数。值得注意的是：一些简单的非线性最小二乘问题通常先作变换将问题转化为线性最小二乘问题求解。

评论：（1）先作变量代换对新变量求最小二乘拟合函数，然后再还原所得近似函数与直接对原变量按最小二乘原则求得拟合函数是不同的。但由于实际计算时，人们主要关心的是问题的简化，就把两者较小的差别忽略了。

（2）以上我们是通过描点观察或经验估计来确定拟合函数的形式，更一般的拟合函数的选择问题，请参考冯康所编著的《数值分析》。 （3）当

n3时，最小二乘法的正则方程组一般是病态的，n愈大病态情形

f(x)与g(x)的内积为

mi1更严重。为了避免求解病态方程组，我们必须引入点集上的正交函数族。

对离散和连续两种情形，通过引进内积与范数的概念，将它们统一起来。在离散情形，我们定义函数

(j,k)ij(xi)k(xi)

在连续情形，则定义函数

f(x)与g(x)的内积为

ba(f,g)(x)f(x)g(x)dx

容易验证以上两种均定义了内积空间。

例5 利用正交函数族求例2所给数据表的最小二乘二次拟合多项式。 解 按式(2.9)和(2.10)计算，得

(x0,0)80(x)1, 1xi(0,0)i0(x1,1)831(x)x, 2xi(1,1)i08(1,1)1xi2(0,0)i0810,i088xi02i0,

3.7510.41667,9i02(x)x20.41667.由式(2.7)，得

8(y,)*0a0yi(0,0)i08(y,)*1a1yixi(1,1)i0818.172312.019149i08.4842x2.26253.75i02i88(y,2)a(2,2)*22y(xii0.41667)i0822(x0.41667)ii00.0378

得最小二乘二次拟合多项式为

**0011*(x)a(x)a(x)a22(x)2.019142.2625x0.0378(x20.41667) 2.00342.2625x0.0378x2性质1（正交性）

0, nm11Pn(x)Pm(x)dx2,nm

2n12例8 设f(x)1x，求[0,1]上的一次最佳平方逼近多项式. 解这是(x)1的情形.取0(x)1,1(x)x, span{1,x}，于是

111(0,0)1dx1,(0,1)xdx,002112(1,1)xdx,03

d0(f,0)得方程组

101x2dx1.147,d1(f,1)x1x2dx0.60901112a01.1471213a0.609 1解出a00.934,a10.426。故

S1*(x)0.9340.426x

平方误差

22f(x),f(x)S1*(x),f(x)(1x2)dx0.426d10.934d00.002601

最大误差

maxf(x)S(x)0x1*max1x20.9340.426x

0x10.066令

h(x)1x0.9340.426x，则

h(x)x221x10.4262，由h(x)0，知x0.426h(0)0.066,h(1)0.054，可得0.066。

用

0.426

结合

{1,x,,x}做基，求最佳平方逼近多项式，当n较大时，Hilbert

n系数矩阵(5.9)是高度病态的，因此直接求解法方程是相当困难的，通常是采用正交多项式做基。

f(x)ex在[1,1]上的三次最佳平方逼近多形式

解 先计算(f(x),Pk(x)) (k0,1,2,3)

例9 求

1(f(x),P0(x))edxe2.35041e

1xx1(f(x),P(x))xedx2e0.7358111321x(f(x),P2(x))xedx0.1431

122153x2(f(x),P3(x))xxedx0.02013

1221于是

*a0(f(x),P0(x))21.1752,*a13(f(x),P1(x))21.1036, *a25(f(x),P2(x))20.3578,a7(f(x),P3(x))20.07046

因此

*S3(x)0.99630.9979x0.5367x20.1761x3

*3均方误差

3(x)2eS(x)1x*323

2*2e2xdxak0.00841k02k1最大误差

3(x)*exS3(x)0.0112

例12 给出有理函数 42x45x3381x21353x1511 R43(x)32x21x157x409用辗转相除法将它化为连分式，并写成紧凑形。

解 用辗转相除可逐步得到

4x264x284R43(x)2x33x21x2157x40942x36(x9)x52x16x7142x36x58 x7x94682x3

x5x7x91kxx展成连分式. 例13 将级数ek0k!利用定理15，得

11xkexkk!k!xk0k01111112341x2x6x24x11x24x436x6112233411xx2x2x6x6x24x1xx4x36x1x1x22x66x241xx2x3xkx1x1x2x3x4xk1

例14考虑函数

ln(1x)的逼近问题。它的Taylor展开式为

kxln(1x)(1)k1,x(1,1] (6.6)

kk1n取部分和

kxTn(x)(1)k1ln(1x)

kk1 另一方面，若对(6.6)式用辗转相除或者定理15可得到ln(1x)的一种连分

式展开

x1x1x2x22xln(1x)12345若取(6.7)式的2,4,6,8项，则可分别得到ln(1x)的以下有理逼近

(6.7)

2x6x3xR11,R222x66xx260x60x211x3R33

236090x36x3x420x630x2260x325x4R44420840x540x2120x36x4若采用同样多项的泰勒展开部分和S2n(x)逼近ln(1x)，并计算ln2的值S2n(x)及Rnn(x)，计算结果见表3.3，R44(1)的精度比S8(1)高出

10万倍，而它们的计算量是相当的，这说明有理逼近好得多，开展某些函数的有理逼近时有必要的。

表3.3 计算结果

2n1 S2n(1) Sln2S2n(1) 1 2 3 4 0.5 0.58 0.617 0.634 0.19 0.11 0.076 0.058 Rnn(1) 0.667 0.692 31 0.693 122 0.693 146 42 Rln2Rnn(1) 0.026 0.000 84 0.000 025 0.000 000 76 f(x)ln(1x)的帕德逼近R(2,2)及R(3,3).

解 由ln(1x)的泰勒展开

121314ln(1x)xxxx

234得c00,c11,c212,c313,c414,。当nm2时，由(6.14)得

例14 求

11b2b123 1b1b121342求得1b1,b216，再由(6.12)得

a00,a11,a212

于是得

12xx26x3x2 R22(x)21266xx1xx6当nm3时，由(6.14)得

111b32b23b141111b3b2b1

345211113b34b25b16331,b2,b3解得b1.代入(6.10)得a00, 252011a11,a21,a3.于是得

601132xxx2360x60x11x60 R33(x)23332136090x36x3x1xxx2520帕德逼近Rnm(x)的误差估计.由(6.12)及(6.14)求得Pn(x),Qm(x)的

,an及b0,b1,,bm，代入则得 数a0,a1,mlnm1f(x)Qm(x)Pn(x)xbkcnm1lkx

l0k0将Qm(x)除上式两端，即得

f(x)Rnm(x)xnm1rlxlQm(x)

l0其中l当

rbkcnm1lk.

k0mx1时可得误差近似表达式

f(x)Rnm(x)r0xnm1，r0bkcnm1k

k0m第4章 数值积分与数值微分

例1 判断求积公式

11f(x)dx9[5f(0.6)8f(0)5f(0.6)]

1的代数精度。

解 记

I(f)f(x)dx111I(f)[5f(0.6)8f(0)5f(0.6)]9因为

1I(1)dx2，I(1)(585)2191I(x)xdx=0111I(x)[50.6805(0.6)]091222I(x)xdx=13122I(x)(50.68050.6)93I(x)x3dx=01312I(x)xdx=15124I(x)(50.36050.36)95414 1I(x3)[5(0.6)305(0.6)3]09I(x)xdx=015151I(x)[5(0.6)505(0.6)5]091266I(x)xdx=1712 633I(x)[5(0.6)05(0.6)]0.24975所以求积公式具有5次代数精度。

例2给定形如

10f(x)dxA0f(0)A1f(1)B0f(0)的求积公

式，试确定系数A0,A1,B0，使公式具有尽可能高的代数精度。

2解 求积公式中有三个参数，因此至少对f(x) 1,x,x应精确成立，即 当f(x)1时，得

当

f(x)x时，得

A0A11dx1

01当

f(x)x21A1B0xdx

021时，得

121A1xdx

03211,A1,B0，于是有 解得A033612110f(x)dx3f(0)3f(1)6f(0)

111333当f(x)x时，xdx，而上式右端为，故公式对f(x)x034不精确成立，其代数精度为2。

例2 根据函数表4-1，用复化梯形公式和复化辛普森公式计算

sinxIdx的近似值，并估计误差。

0x解 将积分区间[0,1]划分成8等分，应用复化梯形公式有

71kI[f(0)f(1)2f()]0.9456909

168k1而如果将积分区间[0,1]划分成4等

1分，应用复化辛普森公式得

341k2k1I[f(0)f(1)2f()4f()] 2448k1k10.9460832与准确值I0.9460831比较，显然用分，应用复化辛普森公式得

341k2kI[f(0)f(1)2f()4f(2448k1k10.9460832

与准确值

I0.9460831比较，显然用复化Simpson公式计算精度较高。

为了利用余项公式估计误差，要求

h

sinxf(x)x的高阶导数，由于

1sinxf(x)cos(xt)dt0x

所以有

1kf(k)1dkk(x)cos(xt)dttcos(xtk0dx02

于是

maxf0x1(k)k(x)maxtcos(xt)dt00x1211ktdt0k11k

由复化梯形误差公式(3.3)得

22h111R8(f)maxf(x)0.00120x11283

由复化辛普森误差公式(3.6)得

bah(4)R4(f)maxf(x)0x1180211160.2711018085 例3 若用复化求积分公式计算积分

1x044Iedx

的近似值，要求计算结果有四位有效数字，

n应取多

大？

解 因为当

0x1时，有

0.3e1ex1

于是

0.3exdx101

要求计算结果有四位有效数字，即要求误差不超过

1410。2又因为

(k)f(x)ex1 x[0,1]2

由式(3.3)得

12111RThf()1041212n2

即

1n10462，开方得

n40.8。因此若用复化梯形公式求积分，

n应等

于41才能达到精度。

若用复化Simpson公

式，由式(3.6)

1h(4)11RSf()180218016n

即。

故

得取

44n1.62n2，即

n2的复

合Simpson公式计算即可达到精度要求，此时区间

[0,1]实际

上应分为4等分，只需计算5个函数值，而复合梯形公式则需计算43个函数值，工作量相差8倍之多。

4龙贝格求积

公式 4.1梯形公式的逐次分半算法

如前所述，复化求积公式的截断误差随着步长的缩小而减少，甚至由给定的

精度可以预先确定步长。但由于误差估计式中出现高阶导数，使问题变得十分困难。实际计算时，我们总是从某个步长出发计算近似值，若精度不满足要求，可将步长逐次分半以提高近似值，直到求得满足精度要求的近似值。

设将区间

[a,b]分为n等分，共有n1个分

点，如果将求积区间再二分一次，则分点增至

2n1个，我们将二分前后两个积分值联系起来加以考虑。注意到每个子区间

[xk,xk1]经过二分只增加了一个分点

xk121(xkxk1)2，用复化梯形公式求得该子区间上的积分值为

hf(xk)2f(xk12)f(xk1)4

注意，这里

bahn代表二分前的步长，将每个子区间上的积分值相加得

hn1hn1T2nf(xk)f(xk1)f(xk124k02k0

即

1hn1T2nTnf(xk12)22k0 (4.1) 这表明，将步长由

h缩小

hT2n为时，2等于Tn的一

半再加新增加节点处的函数值乘以当前步

长。 由复合梯形误差公式

h3n1ITnf(k)，k(xk,xk1)12k0

知

n1ITn11blim2limhf(k)f(h0h12h0k012a1f(b)f(a)12

所以当

n充

2分大时，有

hITnf(b)f(a)121ba12n

及

2f(b)f(a)21baIT2n122n

得

f(b)f(a)IT2n1ITn4

有

1IT2n(T2nTn)3 (4.2)

这说明，若把

T2n作为积分值I的近

似值，误差大约

为

(T2nTn)3。

算法4.1 1．输入

a,b,f(x),

2．置

bam1,h,T0h[f(a)f(b)]23．置

F0，对

k1,2,,2m1

FFf(a(2k1)h) 4

．

1TT0hF25．若

TT03，

输

出

IT，停

机；否则

hm1m,h,TT02，转3。

4.2 李查逊

（Richardson）外推法 假设用某种数值方法求量

I的近似值，

一般地，近似值是步长

h的函数，记为

I1(h)，相应

的误差为

II1(h)1hp12hp2 (4.3)

其

中

khpkpki(i1,2,),0p1p2是与h无关

的常数。若用

h则得

代替

(4.2)中的

h，

k(h)pkkhpkpkII1(h)1(h)p11h (4.4) 式(4.4)减去式(4.3)乘以p1，得

p1p1II1(h)[II1(h)]2(p2p1)hp23(p3p1)hp3

p1取

k(pkp1)hpk

满足1，以

p1除

1上式两端，得

I1(h)p1I1(h)p3p2Ibhbh23p11 (4.5)

其中

pip1p1bkhpkbi2()(1)(i2,3,)仍与h无关。

令

I1(h)p1I1(h)I2(h)1p1

由式(4.5)，以

I2(h)作为I的近似值，

其误差至少为

p2，

O(h)因此I2(h)收敛于I的

速

度

比

I1(h)快。不

断重复以上作法，可以得到一个函数序列

Im1(h)pm1Im1(h)Im(h)，m2,3,pm11 (4.6)

I(h)近似I，误差

以m为

IIm(h)O(h)。随着m的

增大，收敛速度越来越快，这就是Richardson外推法。

4.3 龙贝格求积公式

由前面知道，复化梯形公式的截断误差

为

pmO(h2)。进

一步分析，我们有如下欧拉—麦克劳林（Euler-Maclaurin）公式：

定理4 设

f(x)C[a,b]，则有

IT(h)1h2h

其中系数

24kh2kk(k1,2,)与h无关。

把李查逊

外推法与欧拉—麦克劳林公式相结合，可以得到求积公式的外推算法。特别地，在外推算法式(4.6)中，取

1,pk2k2，

并

记

T0(h)T(h)，则有

4Tm1(h2)Tm1(h)Tm(h),m1,2,m41 (4.7)

经

过

mm(m1,2,)次加速后，余项便取下列形式：

Tm(h)I1h (4.8) 上述处理方法通常称为李查逊（Richardson）外推加速方法。

为研究Romberg求积方法的机器实现，引入记号：

2(m1)2h2(m2)以T0二分

(k)表示

k次后

(k)求得的梯形值，且以Tm表示序列

(k)T0的m次加速值，

则依以上递推公式得到

m(k1)(k)4TT(k)m1m1Tm(h),k0,1,2,m41

称为龙贝格求积算法。

Romberg公式的计算过程见下表4-2

表4-2

k

表4-1

k xk sinxkf(xk)xk k xk sinxkf(xk)xk 0 0 1 1 18 0.9973978 2 14 0.9896158 3 38 0.9767267 5 58 0.9361556 6 34 0.9088516 7 78 0.8771925 8 1 0.8414709 4 120.9588510 0 1 2 T0(k) T1(k) T2(k) T3(k) T4(k) ba T0(0)  ba(1)(0) T0 T1 2 3 4 ba(2)(1) T0 T1 4ba(3)(2) T0 T1 8ba(4)(3) T0 T1 16 1T2(0) T2(1) T2(2) (0)T 3 T(0) (1)T 34

例4 用Romberg算法计算积分Ix32dx。

0解 利用逐次分半算法(4.1)和Romberg算法(4.7)，计算结果见表4-3。

T(0)0T0(1)T(2)0T0(3)1[f(0)f(1)]0.50000021(0)11T0f()0.4267772221(1)0.513 T0[f()f()]0.40701822441(2)0.251357T0[f()f()f()f()]2288880.401812 表4-3

k 0 1 2 3 4 5 T0(k)

0.500000 0.426777 0.407018 0.401812 0.400463 0.400118 1T1(k)

0.402369 0.400432 0.400077 0.400014 0.400002 T2(k)

0.400302 0.400054 0.400009 0.400002 T3(k)

0.400050 0.400009 0.400002 T4(k)

0.400009 0.400002 T5(k)

0.400002 1例5计算积分

0.2x2dx，若用复合辛普森法(3.5)式计算结果见表4-6(此处hn即为公式中的h，积分精确值为4）

表4-6 计算结果 n hn 1 2 3 4 5

0.8 0.4 0.2 0.1 0.05

Sn 4.948148 4.187037 4.024218 4.002164 4.000154

SnSn1 0.761111 0.162819 0.022054 0.002010

1 1计算到SnSn10.02为止，此时I(f)0.2x2dx的近似值S5[0.2,1]4.000154，若再用龙贝格法则得到

S5S4RS[0.2,1]S54.00002

15整个计算是将[0.2,1]做32等分，即需要计算33个f(x)的值。现在若用自

适

应

积

分

法

，

当

h20.4时有

S2[0.2,0.6]3.51851852,S2[0.6,1]0.66851852，由于

S2S2[0.2,1]S2[0.2,0.6]S2[0.6,1]4.187037,

S1S20.761111

h2大于允许误差0.02，故要对[0.2,0.6]及[0.6,1]两区间再用h3做积

2分。先计算的积分[0.6,1]S3[0.6,0.8]0.41678477,S3[0.8,1]0.25002572。

由于

S2[0.6,1](S3[0.6,0.8]S3[0.8,1])0.668518520.666810490.001708小于允许误差0.01，故在[0.6,1]区间的积分值为

1RS[0.6,1]0.66681049(0.666810490.66851852) 150.66669662下面再计算子区间[0.2,0.6]的积分，其中

h2S2[0.2,0.6]3.51851852，而对h3可求得

2S3[0.2,0.4]2.52314815,S3[0.4,0.6]0.83425926

由于

S2[0.2,0.6](S3[0.2,0.4]S3[0.4,0.6])0.161111

大于允许误差0.01，因此还要分别计算[0.2,0.4]及[0.4,0.6]的积分，当

h3h4时可求得

2S4[0.4,0.5]0.50005144,S4[0.5,0.6]0.33334864

而

S3[0.4,0.6](S4[0.4,0.5]S4[0.5,0.6])0.000859

小于允许误差0.005，故可得[0.4,0.6]的积分近似

RS[0.4,0.6]0.8333428

而对区间[0.2,0.4]，其误差S3[0.2,0.4]S4[0.2,0.4]不小于0.005，

故还要分别计算

[0.3,0.4]的积分，其中

h4S4[0.3,0.4]0.83356954，当h5可求得

2S5[0.3,0.35]0.47620166,S5[0.35,0.4]0.35714758

及

[0.2,0.3]且

S4[0.3,0.4](S5[0.3,0.35]S5[0.35,0.4])0.000220

小于允许误差0.0025，故有

RS[0.3,0.4]0.83333492

最后子区间[0.2,0.3]的积分可检验出它的误差小于0.0025，且可得

RS[0.2,0.3]1.666686

将以上各区间的积分近似值相加可得

I(f)RS[0.2,0.3]RS[0.3,0.4]RS[0.4,0.6]RS[0.6,1]4.00005957

它一共只需计算17个例9 确定求积公式

1f(x)的值。

的系数A0,0xf(x)dxA0f(x0)A1f(x1)

A1及节点x0,x1，使它具有最高代数精度。

解 具有最高代数精度的求积公式是高斯型求积公式，其节点为关于权函数

(x)x的正交多项式零点x0及x1，设(x)(xx0)(xx1) x2bxc，由正交性知(x)与1及x带权正交，即得

于是得

10x(x)dx0,10xx(x)dx0

222222bc0及bc0 753975105,c，即 由此解得b921105(x)xx

921令(x)0，则得x00.289949,x10.821162

2由于两个节点的高斯型求积公式具有3次代数精确度，故公式对

f(x)1,x，精确成立，即

当f(x)1时

A0A1当

10f(x)x时

12xdx

32A0x0A1x1xxdx

05由此解出A00.277556,A10.389111

,n)上下面讨论高斯求积公式的余项。设在节点xk(k0,1,f(x)的2n1次Hermite插值多项式为H(x)，即

n1(xk)f(xk),k0,1,,n H2n1(xk)f(xk),H2由Hermite余项公式

f(2n2()2f(x)H(x)n1(x)

(2n2)!有

R(f)(x)f(x)dxAkf(xk)ak0nbn(x)f(x)dxAkH(xk)ak0bb(x)f(x)dx(x)H(x)dx

aab(x)f(x)H(x)dxabf(2n2()2(x)n1(x)dxa(2n2)!例10 用4点(n3)的高斯-勒让德求积公式计算

b解 先将区间[0,2Ixcosxdx

202]化为[1,1]，由（6.13）式有

1I()(1t)cos(1t)dt

144根据表4-6中n3的节点及系数值可求得

32IAkf(xk)0.467402 (准确值I0.467401例11 用5点(nk03)

5)的高斯-切比雪夫求积公式计算积分

1exIdx

211xx(2n)(x)ex，当n5时由公式(6.14)可得 解 这里f(x)e，fIe5k15cos2k1103.977463

9由余项(6.16)式可估计误差

R[f]210!9e4.610

例12 用高斯-拉盖尔求积公式计算

的近似值。 解 取n0exsinxdx

1，查表得

x00.58578644,x13.41421356.A00.85355339,A10.14644661故

，

若取n0exsinxdxA0sinx0A1sinx10.43246

02，可得exsinxdx0.49603；

若取n5，可得exsinxdx0.50005。

0而准确值确。

0exsinxdx0.5，它表明取n5的求积公式已相当精

区间为(,)，权函数

(x)e2x2的正交多项式为埃尔米特多项式

ndnxx2Hn(x)(1)ee,n0,1,ndx,

对应的高斯型求积公式

exf(x)dxAkf(xk)

k02n称为高斯-埃尔米特求积公式。节点xk 特多项式的零点，求积系数为

(k0,1,,n)为n1次埃尔米

Ak2(n1)!n11(xk)][Hn2高斯-埃尔米特求积公式的节点和系数可见表4-9

表4-9 高斯-埃尔米特求积公式的节点和系数 n 0 1 2 3 4 xk 0 ±0.707106781 ±1.2247448710 ±1.650680124 ±0.524647623 ±2.020182871 ±0.9585724650 ±2.350604974 ±1.335849074 ±0.436077412 ±2.651961357 ±1.673551629 ±0.8162878830 Ak 1.772453851 0.886226926 0.295408975 1.181635901 0.081312835 0.804914090 0.019953242 0.393619323 0.945308721 0.004530010 0.157067320 0.724629595 0.0009717812 0.0545155828 0.425607253 0.810264618 5 6 公式(6.20)的余项为

(n1)!(2n2)R[f]n1f(),(,) (6.22)

2(2n2)!例13 用两个节点的高斯-埃尔米特求积公式(6.20)计算积分解 先求节点

ex2x2dx

x0,x1，由

H2(x)4x22，其零点为

22x0,x1，由(6.21）式可求得

22A0A1于是

x222

2222exdx2[(2)(2)]2

2高斯型求积公式代数精度为3，故对f(x)x求积公式精确成立，从而得

2ex2xdx22

例14 用复合辛普森公式求二重积分

的近似值。 解 取N1.51.41.0ln(x2y)dydx

2,M1，即h0.3,k0.5得

21.51.41.0ln(x2y)dydx22k2[ln(x2)dx4ln(x2.5)dxln(x3)dx]1.41.461.40.50.3[ln3.44(ln3.55ln3.85)2ln3.7ln4]660.51.2[ln3.94(ln4.05ln4.35)2ln4.2ln4.5]660.50.3[ln4.44(ln4.55ln4.85)2ln4.7ln5]660.42955244

此积分的真值是0.4295545265（保留小数后10位） 对二重积分(7.1)式也可用其他求积公式计算，特别是为了减小函数值计算可采用高斯求积公式。

2的高斯求积公式求例14中的二重积分

解 先将区域R{(x,y)|1.4x2,1.0y1.5}变换为区域R{(u,v)|1u,v1}，其中

11u(2x3.4),v(2y2.5)

0.60.5或等价于

例15 用nx0.3u1.7,y0.25v1.25

于是有

I对于u,v取n2.01.41.51.0ln(x2y)dydx1110.0751

ln(0.3u0.5v4.2)dvdu2时的高斯求积公式节点及系数，即

u0v00.774596662,u1v1058

u2v20.774596662,A0A2,A199用n2的高斯求积公式计算积分I可得

I0.075111221ln(0.3u0.5v4.2)dvdu0.075AiAjln(0.3ui0.5vj4.2)

i0j00.42955453这里只需计算9个函数值。而例14中需求15个函数值，这里的精度也比例14高，达到8位有效数字。

对于非矩形区域的二重积分，只要化为累次积分，也可类似矩形域情形求得其近似值，如二重积分

I用辛普森公式可转化为

bbad(x)c(x)f(x,y)dydx

k(x)I[f(x,c(x))4f(x,c(x)k(x))f(x,d(x))]dx

a3d(x)c(x)其中k(x)，然后再对每个积分使用辛普森公式，则可求得

2积分I的近似值。

例0.075如，用中点公式求f(x)x在x2处的一阶导数

2h2h G(h)2h设取4位数字计算，结果见表4-10（导数的准确值f(2)0.353553）

表4-10 计算结果

h G(h) 1 0.5 0.1 0.3660 0.3564 0.3535 h 0.05 0.01 0.005 G(h) 0.3530 0.3500 0.3500 h 0.001 0.0005 0.0001 G(h) 0.3500 0.3000 0.3000 从表4-10中看到h0.1的逼近效果最好，如果进一步缩小步长，则逼近效

果反而越差。这是因为当f(ah)及f(ah)分别有舍入误差1和2。若令max{1,2}，则计算f(a)的舍入误差上界为

12(f(a))f(a)G(a)

2hh表明h越小，舍入误差

(f(a))越大，故它是病态的。用中点公式(8.1)计算

f(a)的误差上界为

h2E(h)M

6h要使误差E(h)最小，步长h应使E(h)0，由

hE(h)M20

3h33可得h3/M，如果h3/M，有E(h)0；如果h33/M，有E(h)0。由此得出h33/M时E(h)最

35/2小。当, f(x)x时，f(x)x8135/2Mmaxx0.07536。假定104，则

1.9x2.18241.510h30.125。与表4-10基本相符。

0.075362x例16 用外推法计算f(x)xe在x0.5的导数。

11(h)(h)111[(h)2e2(h)2e2]，当解 令G(h)2h22h0.1,0.05,0.025时，由外推法表4-11可算得

G(0.1)0.4516049081,G(0.05)0.4540761693G1(h)0.4548999231hG(0.025)0.4546926288G1()0.4548981152G20.4548979942

f(0.5)的精确值0.454897994，可见当h0.025时用中点微分公式

只有3位有效数字，外推一次达到5位有效数字，外推两次达到9位有效数字。

第五章 线性代数方程组的数值解法

例1 解方程组

102x1x21 x1x22(1) 求精确解； (2) 在

10，n2的浮点机上用Gauss消去法求解。

11,x21。 解 (1) 容易求得方程组的精确解为x119999(2) 在所给浮点机上原方程组为

1110.1010x0.1010x0.101012 1110.1010x10.1010x20.2010由于

m210.10101(0.10101)0.10103

消去第二式的x1，得

(0.101010.10103)x20.201010.10103

对阶，出现大数“吃掉”小数，结果有

33 21解得2，回代第一式得1。

与精确解相比较，已无精确度可言。产生不准确的原因是主元素太小，致使

0.1010x0.1010x0x0.1010m21很大。改变上述状况的办法是将方程组的一、二两式对换，得

0.1010x10.1010x20.2010 1110.1010x10.1010x20.1010然后再使用Gauss消去法，此时

12111111m0.1010(0.1010)0.1010

(0.10100.1010)x20.101010.20101

11于是得近似解x10.1010,x20.1010。

此例表明，高斯消去法解方程组时，小主元可能带来严重的后果，因此应尽量避免小主元的出现；另一方面，通过对换方程组中方程的次序或改动变量次序，选择绝对值大的元素做主元，可减少舍入误差，提高计算精度。

11

例2 求矩阵

223A477

245的三角分解。

解 按式(3.2)

u11a112,u12a122,u13a133l21a21u11422,l31a31u11221u22a22l21u127223u23a23l21u137231l32(a32l31u12)u22[4(1)2]32u33a33(l31u13l32u23)5[(1)321]6所以

100223A210031

121006紧凑格式：

(a11) u11 (a12) u12 (a13) u13 (a21) l21 (a22) u22 (a23) u23

… … … …

(a1n) u1n (a2n) u2n

(a31) l31 (a32) l32 (a33) u33

(a3n) u3n

(an1) ln1 (an2) ln2 (an3) ln3

例2中矩阵

(ann) unn

A的三角分解按紧凑格式计算，结果见下表

(2) 2

(3) 3

(2) 2

4(7) 7223 (7) 7231 (4) 2

212(2) 1 (4) [4(1)2]2 (5) 521(1)36 32Axb的系数矩阵已进行三角分解，ALU，则解

方程组Axb等价于求解两个三角形方程组Lyb,Uxy。

第一步：先求解下三角方程组Lyb

如果线性方程组

1lLy21ln1得解

1lnn1y1b1 1ynbny1b1k1 ykbklkjyj, k2,3,,nj1第二步：再求解上三角方程组Uxy

u1nx1y1u11u12u22u2nx2y2  unnxnyn得解

xnynunnn1x(yux), kn1,n2,kkkjjukkjk1例3 用选主元的直接三角分解法解方程组Axb，其中

,1

113A246，b(1,4,1)T,x(x1,x2,x3)T

492解(1) 对A作选主元的LU分解，中间结果冲掉A相应位置元素，数组

Ip(3)记录主行号。

第一步分解：a11S11,a21S22，a31S34，因为S3maxSi，于是有

1i3492492r1r3分解A2461246

113141315第二步分解：a22S2,a23S3，因

24S3maxSi，故有

2i3为

41r2r3A4124921分解1341462492第三步155分解 442124259542525 26于是

1L141492，U5452

122514由选主元过程，Ip(1)3,Ip(2)3,Ip(3)3。 001PI23I13100，Pb(1,1,4)T010 (2) 解LyPb，得y(1,34,3210)T；解Uxy，x(0.8,1,2.4)T

例6 用改进平方根法求解方程组 解 容易验证，系数矩阵

121A250

1014为对称正定阵。按式(4.6)计算分解式，得

d1a111t21a212,t31a311,l21t21d12,l31t31d11 d2a22t21l215221

t32a32t31l210122l32t32d22 d3a33t31l31t32l321411(2)(2)9按式(5.7)计算方程组的解，得

得

y1b14,y2b2l21y17241y3b3l31y1l32y21514(2)(1)9 x3y3d31,x2y2d2l32x31(2)11x1y1d1l21x2l31x342111所以方程组的解为1

x1,x21,x31。

例6 容易看出，方程组

x1x22 x11.00001x22的解为x12,x20。而方程组

x1x22 x11.00001x22.00001的解为x11,x21。

例7 计算例6方程组系数矩阵的条件数。 解 系数矩阵为

11A

51110其逆矩阵为

51101A51010 510

5于是有

cond(A)AA14105(2105)(21051)条件数很大，所以方程组是“病态”的。

例8 已知希尔伯特（Hilbert）矩阵

11Hn21n计算

12131n11n1n1

12n1H3的条件数。

1解 Hn的逆矩阵Hnaij的元素是

nn(1)ij(ni1)!(nj1)!(ij1)(i1)!(j1)!(ni)!(nj)!

2aij所以

, 1i,jn3630112139H3121314，H3136192180

13141530180180（1）计算H3的件数cond(H3)。

H3所以

116,H31408

cond(H3)748

同样，可算得

cond(H6)2.610当

7，cond(H7)9.85108

n愈大时，Hn病态愈严重。

（2）考虑方程组

TH3x(116,1312,4760)b

的扰动方程（H3及b的元素取3位有效数字）

1.000.5000.333x1x11.830.5000.3330.250xx1.08

220.3330.2500.20xx0.78333简记为

(H3H3)(xx)bb

方程

H3xb与它的扰动方程的解分别为：

x(1,1,1)T，

xx(1.089512538,0.487967062,1.491002798)T

于是

x(0.0895,0.5120,0.4910)而

T，

xx51.2%

00.0003300.00333，b0.00333

H300.0003300.000330.0003300有

H3所以

0.00033，b0.00333

H3H3b0.003330.000330.02%,0.182%116b116

这表明

H3与b的相对误差不超过0.3%，而引起解的相对误差超过

450%。例9 设

4x110110x 1122简记为Axb，计算cond(A)。

解

4110411101A,A4

1011111411041A110，A4

101有

cond(A)AA1（110）410

410142当用列主元消去法求解时（计算到三位有效数字）

1104104 (A|b)4401010于是得到很坏的结果：

x21,x10。

41100111取D，考察等价方程组DAxDb的系数矩101阵ADA的条件数。

11041111A ,A441101101144，大大改善了系数矩阵的条件数。 则cond(A)4110再用列主元消去法求解

Ax=b，得到

10411112112(A|b) 4  

1101111210从而得到较好的计算解：

x21,x11。

例1 用直接三角分解法解方程组

7.0006.990x134.974.0004.000 x20.00 2解 因方程组的系数矩阵的两行（列）几乎线性相关，所以方程组可能是病态的。事实上

cond(A)A1方程组的确是病态的。 利用迭代改善法： （1）实现A=A275.1413.993849

LU分解得

07.0006.9901A  LU

0.0060.571410（2）解Lyb及Uxy得

y(34.97,0.02)T,x(1.667,3.333)Tx(0)

(0)(0)TrbAx(0.00333,0)（3）计算（用双倍字长计算）

(0)（4）求解LUzr(0)得

z(0)(0.3171,0.3172)T

（5）计算改善解

（6）计算

x(1)x(0)z(0)(1.986,2.999)T

x(1)x(0)0.317。重复（3）—（6）步做法，得

r(1)bAx(1)(0.0002,0)T (1)Tz(0.00241,0.01667) x(2)x(1)z(1)(1.986,2.999)T

1(2)(1)xx0.017101

21上述迭代改善过程还可以继续下去，如果给定101，则x(2)就是满

2*T足精度要求的近似解。方程组的准确解x(2,3)。