计算方法作业-1(2013研究生)

《计算方法》课程 No.1

专业名称： 所在学院： 导师姓名： 姓 名： 学 号：

完成时间： 作业成绩： 任课教师签字：

第一次作业：

请同学们完成下列作业（提交时间：2013.10.25）

1. 已知x1 1.03 0.01，x2 0.45 0.01，试求用公式y x1 生的绝对误差和相对误差。

解 e(x1) 0.01，e(x2) 0.01，

2

1x2

e计算y时产2

y1x2 y

2x1， e，故绝对误差为

x22 x1

e(y) d(y)

y y

e(x1) e(x2) x1 x2

2.06 0.01 0.5 e0.45 0.01 0.0284416

相对误差为er(y)

2. 建立积分In 解 由

e(y)dy0.0284416

0.015415 20.45yy1.03 0.5 e

10

xn

dxn 0,1, ,20的递推关系式，并研究它的误差传递。 x 5

In 5In 1

和

1xn 5xn 11

dx xn 1dx 005 xn1

I0

可建立下列递推公式

1

dx ln6 ln50x 5

1

I0 ln6 ln5 1I 5In 1n n

(n 1,2, ,20)

(*)

计算出I0后，由递推关系式可逐次求出I1,I2, ,I20的值。但在计算I0时有舍入误差，因

*

此在使用递推关系式中，实际算出的都是近似值In(n 1,2, ,20)。即

I0* I0 e0

*1* In 5In 1

n

误差是如何传递的？

(n 1,2, ,20)

设I0有误差e0，假设计算过程中不产生新的舍入误差，则由(*)式可得

**

en In In 5In 1 5In(n 1,2, ) 1 5en 1

*

从而有

en ( 5)ne0

即原始数据I0的误差e0对第n步的影响使该误差扩大了5倍。当n较大时，误差将淹没真值，因此递推公式(*)是数值不稳定的。

现在从另一方向使用这一公式

*

n

1 1

In 1 In (n 20,19, ,1) (**)

5 n

****

只要给出I20的一个近似值I20，即可递推得到I19，类似于上面的推导可得 ,I18, ,I0

en 1

每递推一步误差缩小到原值的

由于x [0,1]时，

1 1

en,e0 en

5 5

n

1

，所以递推公式(**)是数值稳定的。 5

xnxnxn

65 x5

所以有估计式

11

In

6(n 1)5(n 1)

于是

11

I20 6 215 21

取

11

I20 0.0087301587

2

I20 0.0087301587

可得另一算法： 11

In 1 ( In)(n 20,19, ,1) 5n

这是数值稳定性好的算法。

3. 用秦九韶法计算P5(x) 6x5 5x4 4x3 3x2 2x 1在x 5处的值，并验证之。 答案：22461

解 Mathematica程序：

Clear[a,S,k,x] a[k_]:=k+1; S[5]=a[5];

S[k_]:=x*S[k+1]+a[k]; x=5; S[0]

运行结果：22461

秦九韶算法也可改用下面的简单竖式算法：

6 5 4 3 2 1

30 175 895 4490 22460

6 35 179 898 4492

用Mathematica语句验证： a[k_]:=k+1

P[x_]:=Sum[a[k]*x^k,{k,0,5}] P[5]

结果：22461

4. 已知函数y 2x的函数表

0.3

（1）试以x2 0，x3 1为节点，建立线性插值多项式L1(x)，用以计算2值，并估计截断误差。

的近似

0.3

（2）以x1 1，x2 0，x3 1为节点，建立抛物线插值多项式L2(x)，用以计算2的近似值，并估计截断误差。

解（1）线性插值多项式为

L1(x)

则

x 1x 0

1 2 x 1 0 11 0

20.3 L1(0.3) 1.3

其截断误差由(2-13)得

R1(x)

M2

(x x2)(x x3) 2

其中M2 maxf"(x)，因为f"(x) (ln2)22x，可得

x2 x x3

maxf"(x) 2(ln2)2 0.9069

0 x 1

于是

R1(0.3)

（2）抛物插值多项式为

0.9069

0.3(0.3 1) 0.09522 2

L2(x)

则

(x 0)(x 1)(x 1)(x 1)(x 0)(x 1)

0.5 1 2 0.25x2 0.75x 1

( 1 0)( 1 1)(0 1)(0 1)(1 0)(1 1)

20.3 L2(0.3) 1.248

其截断误差为

R2(x)

M3

(x x1)(x x2)(x x3) 6

3

其中M3 maxf'''(x) 2(ln2) 0.6660，

1 x 1

故有

R2(0.3)

5. 已知函数表

0.6660

(0.3 1)(0.3 0)(0.3 1) 0.03030 6

求出Lagrange插值多项式，并计算x 1.2处的y的近似值。

解 Mathematica程序如下：

r0[x_,x0_,a_,b_,c_]:=(x-a)(x-b)(x-c)/((x0-a)(x0-b)(x0-c))

L3[x_]:=r0[x,1,2,4,5]*16+r0[x,2,1,4,5]*12+r0[x,4,1,2,5]*8+r0[x,5,1,2,4]*9 L3[x]//N Simplify[%] L3[1.2]//N 插值多项式为

-1.33333 (-5. + x) (-4. + x) (-2. + x) + 2. (-5. + x) (-4. + x) (-1. + x) -1.33333 (-5. + x) (-2. + x) (-1. + x) + 0.75 (-4. + x) (-2. + x) (-1. + x) 插值多项式的简式为

23

74.8333x3 0.08333x33 0.08333x33 20.666

插值为 15.1307

(xi－x)2li(x)=0，其中li(x)(i=0,1, ,5)为插值基函数。 6. 证明∑

i=0

5

5

5

证明

5

(xi－x)∑

i=0

2

li(x)=∑(xi2－2xxi+x2)li(x)

i=0

5

5

xili(x)－2x∑xili(x)+x ＝∑

2

i=0

i=0

2

li(x)=x2－2xx+x2＝0 ∑

i=0

7. 设f(x) C4[a,b]，已知插值节点 a<x0< <xn<b且xi－xi－1＝h,i=1,2,3，

(b－a)2

maxf''(x) 证明：（1）maxf(x)－L1(x)≤

a≤x≤ba≤x≤b8

（2）maxf(x)－L2(x)≤

x0≤x≤x2

3

hmaxf'''(x) x0≤x≤x227

14

hmaxf(4)(x) （3）maxf(x)－L3(x)≤x0≤x≤x324x0≤x≤x2

证明 （1）由插值余项得

R1(x)=f(x)－L1(x)=

f''(ξ)

(x－a)(x－b) a<ξ<b 2!

又 max(x－a)(x－b)=(x－a)(x－b)x=a+b

a≤x≤b

2

(b－a)2

=

4

(b－a)21

maf''(x) 故 maxf(x)－L1(x)≤max(x－a)(x－b)maf''(x)≤

a≤x≤ba≤x≤ba≤x≤b2a≤x≤b8

当f(x)=x时，等式成立。

（2）由插值余项得

R2(x)=f(x)－L2(x)=

f'''(ξ)

(x－x0)(x－x1)(x－x2) x0<ξ<x2 3!

令 y=x－x0,则有 y－h=x－x1，y－2h=x－x2，0≤y≤2h， 于是

x0≤x≤x2

max

f(x)－L2(x)≤

1

maxy(y－h)(y－2h)maxf'''(x)

x0≤x≤x26x0≤x≤x2

而 g(y)=y(y－h)(y－2h)在0≤y≤2h上的最大值为

从而有

x0≤x≤x2

23

h 9

max

f(x)－L2(x)≤

＝

1

maxy(y－h)(y－2h)maxf'''(x)

x0≤x≤x26x0≤x≤x2

233

h9

x0≤x≤x2

maxf'''(x)

（3）由插值余项得 R3(x)=f(x)－L3(x)=

f

(4)

(ξ)

(x－x0)(x－x1)(x－x2)(x－x3) x0<ξ<x3 4!

令 y=x－x0,则有 y－h=x－x1，y－2h=x－x2，y－3h=x－x3，，0≤y≤3h， 于是

x0≤x≤x3

max

f(x)－L3(x)≤

1

maxy(y－h)(y－2h)(y－3h)maxf(4)(x)

x0≤x≤x34!x0≤x≤x3

==

1

max(y2－3hy)(y2－3hy+2h2)maxf(4)(x)

x0≤x≤x34!x0≤x≤x3

1

4!

9

－h2≤u≤04

max

u(u+2h2)max

x0≤x≤x3

f(4)(x) 其中u y2 3hy

=

14

hmaxf(4)(x) 4!x0≤x≤x3

8. 已知19个点的坐标如下表，根据表中数据求一插值多项式，并由此计算x 11.75点的函数值的近似值。

数据表：

10 2.30259 45

15 2.70805 50

20 2.99573 55

25 3.21888 60

30 3.4012 65

35 3.55535 70

40 3.6888 75

3.80666 80 4.38203

3.91202 85 4.44265

4.00733 90 4.49981

4.09434 95 4.55388

4.17439 100 4.60517

4.2485

4.31749

解

已知数据表：

10 2.30259 45 3.80666 80 4.38203

15 2.70805 50 3.91202 85 4.44265

20 2.99573 55 4.00733 90 4.49981

25 3.21888 60 4.09434 95 4.55388

30 3.4012 65 4.17439 100 4.60517

35 3.55535 70 4.2485

40 3.6888 75 4.31749

Mathematica程序：

A=Table[{x,Log[x]},{x,10,100,5}]//N

g1=ListPlot[Table[A],Prolog->AbsolutePointSize[6], DisplayFunction->Identity]; Interpolation[A,InterpolationOrder->1]

g2=Plot[%[x],{x,10,100},PlotStyle->Thickness[0.05],DisplayFunction->Identity] Show[g1,g2,DisplayFunction->$DisplayFunction] %%%[11.75]

运行结果：

2.

所求插值f(11.75) 2.445。

)的近9. 已知函数y cosx的函数表如下，试用4次等距节点插值公式计算f(0.048

似值，并估计误差。

解 先根据表2-3构造差分表。

用Newton向前插值公式计算f(0.048)的近似值，取x 0.048，h 0.1，

x 0

0.48，用差分表的上半部分差分，得 h

0.48 (0.48 1)

N4(0.048) 1.0000 0.48 ( 0.00500) ( 0.00993)

2

0.48 (0.48 1)(0.48 2) 0.00013

3!

t

0.48 (0.48 1)(0.48 2)(0.48 3)

0.00012

4!

0.99885 cos0.048

误差估计由式(2-38)可得

R4(0.048)

M5

t(t 1)(t 2)(t 3)(t 4)h5 1.5487 10 7 5!

其中M5 sin0.6 0.565。 用Mathematica程序计算如下 x[k]=k*0.1;

y[0]=1.0;y[1]=0.995;y[2]=0.98007;y[3]=0.95534; y[4]=0.92106;y[5]=0.87758;y[6]=0.82534; f1[i_]:=y[i+1]-y[i] f2[i_]:=f1[i+1]-f1[i] f3[i_]:=f2[i+1]-f2[i] f4[i_]:=f3[i+1]-f3[i] f5[i_]:=f4[i+1]-f4[i] f6[i_]:=f5[i+1]-f5[i]

A={{y[0],y[1],y[2],y[3],y[4],y[5],y[6]},{0,f1[0],f1[1],f1[2],f1[3],f1[4],f1[5]}, {0,0,f2[0],f2[1],f2[2],f2[3],f2[4]},{0,0,0,f3[0],f3[1],f3[2],f3[3]}, {0,0,0,0,f4[0],f4[1],f4[2]},{0,0,0,0,0,f5[0],f5[1]},{0,0,0,0,0,0,f6[0]}}; Transpose[A]//N; N[MatrixForm[%],4]

a[0]=y[0];a[1]=f1[0];a[2]=f2[0];a[3]=f3[0];a[4]=f4[0];a[5]=f5[0];a[6]=f6[0];

NN[x]=Sum[a[k]*Product[(t-j),{j,0,k-1}],{k,0,6}]//N N[Expand[%],3]//Chop %/.t->0.048

运行结果：

1. 0 0 0 0 0 0 0.995 -0.005 0 0 0 0 0 0.9801 -0.01493 -0.00993 0 0 0 0 0.9553 -0.02473 -0.0098 0.00013 0 0 0 0.9211 -0.03428 -0.00955 0.00025 0.00012 0 0 0.8776 -0.04348 -0.0092 0.00035 0.0001 -0.00002 0 0.8253 -0.05224 -0.00876 0.00044 0.00009 -0.00001 0.00001

- 0.005 t - 0.00993 (-1. + t) t + 0.00013 (-2. + t) (-1. + t) t

+ 0.00012 (-3. + t) (-2. + t) (-1. + t) t - 0.00002 (-4. + t) (-3. + t) (-2. + t) (-1. + t) t + 0.00001 (-5. + t) (-4. + t) (-3. + t) (-2. + t) (-1. + t) t

化简得： 1 0.00279t 0.00526t2 0.00354t3 0.00117t4 0.0000t6

1.00012 0.99942

) 1.00012。 即 f(0.048

综合实习题（阅读）

（1991年美国大学生数学建模竞赛A题）

某些州的用水管理机构需估计公众的用水速度（单位：G/h）和每天总用水量。许多供

水单位由于没有测量流入或流出水箱流量的设备，而只能测量水箱中的水位（误差不超过5％）。当水箱水位低于最低水位L时，水泵开始工作将水灌入水箱，直至水位达到最高水位H为止。但是依然无法测量水泵灌水流量，因此，在水泵工作时无法立即将水箱中的水位和水量联系起来。水泵一天灌水1~2次，每次约两小时。试估计在任一时刻（包括水泵灌水期间）t流出水箱的流量f(t)，并估计一天的总用水量。

【注1】问题中使用的长度单位为E（英尺，1E＝30.24cm），容积单位是G（加仑，1G＝3.785L）。

【注2】水箱是直径为57E，高为40E的正圆柱体。当水位落入27E以下，水泵自动启动把水灌入水箱，当水位回升至35.5E时，水泵停止工作。

某镇中某一天的真实用水数据见表1。

tt0={0, 3316, 6635, 10619, 13937, 17921, 21240, 25202, 28543, 32284, 35932, 39332, 39435,

43318, 46636, 49953, 53936, 57254, 60574, 64554, 68535, 71854, 75021, 79254, 82649, 85968, 89953, 93270}; MatrixForm[%];

vv0={3175, 3110, 3054, 2994, 2947, 2892, 2850, 2795, 2752, 2697, 0, 0, 3550, 3445, 3350,

3260, 3167, 3087, 3012, 2927, 2842, 2767, 2697, 0, 0, 3475, 3397, 3340}; MatrixForm[%]

LL0=Table[{tt0[[k]],vv0[[k]]},{k,1,28}];

表1 假设

（1） 影响水箱流量的惟一因素是该区公众对水的普通要求，表1中数据反映该镇在通

常情况下一天的用水量（不包括任何非常情况，水位高度、大气情况、温度变化等物理因素对水的流速均无影响）。

（2） 水泵的灌水速度为常数。

（3） 中水箱中流出水的最大流速小于水泵的灌水速度（不让水箱中的水用尽）。 （4） 公众对水的消耗量以全天的活动为基础，可以认为每天的用水量分布均相似。 （5） 水箱的水流速度可用光滑曲线来近似。

引入记号如下：V－水的容积

Vi－时刻ti（h）水的容积 f(t)－流出水箱的水的流速（

G/h） p－水泵的灌水速度（G/h） 为讨论方便，将表1中数据的单位变换为：

时间单位：h

水量： V r2h，单位：103G

变换结果见表2

表2 为求函数f(t)，先求相邻时间区间的中点i 及在时间区间内水箱中流出的水的平均速度 i

tt2=Table[(tt1[[i+1]]+tt1[[i]])/2,{i,1,27}]; MatrixForm[%];

vv2=Table[(vv1[[i]]-vv1[[i+1]])/(tt1[[i+1]]-tt1[[i]]),{i,1,27}]*10^3; MatrixForm[%];

LL2=Table[{tt2[[k]],vv2[[k]]},{k,1,27}];

结果见表3

ti 1 ti

2vi vi 1ti 1 ti

i 1,2,..27 i 1,2,..27.

表3

作出时间t与平均水流量v的散点图

g1=ListPlot[LL2,Prolog->AbsolutePointSize[10]];

求出拟合曲线并观察拟合效果 ft=Fit[LL2,Table[t^k,{k,0,8}],t] g2=Plot[ft,{t,0,26}] Show[g1,g2]

14217.2 - 1881.35 t - 144.448 t2 + 217.299 t3- 54.5544t4 + 6.38979t5 - 0.379256 t6 + 0.0109854 t7 - 0.000123258 t8

由图可见，将t1 0.460556小时和t2 24.460556小时代入f(t)得

f(t1) 13339，f(t2) 10933.3 f(t2) f(t1) 2395.3

一天内用水总量近似等于

Q

24.460556

0.460556

f(t)dt＝275006

Integrate[ft,{t,0.460556,24.460556}]//N

若按常规每千人的用水量为100000G/h计算，因此估计该地区大约有2780人。

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