《计算方法》课程 No.1
专业名称： 所在学院： 导师姓名： 姓 名： 学 号：
完成时间： 作业成绩： 任课教师签字：
第一次作业：
请同学们完成下列作业（提交时间：2013.10.25）
1. 已知01.003.11±=x ，01.045.02±=x ，试求用公式2
2
12
1x e x y +=计算y 时产生的绝对误差和相对误差。
解 01.0)(1=x e ，01.0)(2=x e ，
11
2x x y
=??，2212x e x y =??，故绝对误差为
0284416
.001.05.001.006.2)()()()(45.022
11=??+?=??+??=
≈e x e x y
x e x y y d y e
相对误差为015415.05.003.10284416
.0)()(45.02=?+=≈=e
y dy y y e y e r
2. 建立积分20,,1,05
10
=+=?
n dx x x I n
n 的递推关系式，并研究它的误差传递。 解 由
n x d x x d x x x I I n n n n n 1
5551011
011
==++=+??--- 和
5ln 6ln 5
1
1
00-=+=?
dx x I
可建立下列递推公式
??
???=-=-=-)
20,,2,1(515ln 6ln 10 n I n I I n n (*)
计算出0I 后，由递推关系式可逐次求出2021,,,I I I 的值。但在计算0I 时有舍入误差，因
此在使用递推关系式中，实际算出的都是近似值)20,,2,1(*
=n I n 。即
??
?
??=-=-=-)
20,,2,1(511
**00*0 n I n I e I I n n
误差是如何传递的？
设*
0I 有误差0e ，假设计算过程中不产生新的舍入误差，则由(*)式可得
),2,1(5551
*
11* =-=+-=-=---n e I I I I e n n n n n n 从而有
0)5(e e n n -=
即原始数据*
0I 的误差0e 对第n 步的影响使该误差扩大了n
5倍。当n 较大时，误差将淹没真值，因此递推公式(*)是数值不稳定的。
现在从另一方向使用这一公式
)1,,19,20(1511 =??
?
??-=-n I n I n n (**)
只要给出20I 的一个近似值*
20I ，即可递推得到*0
*18*19,,,I I I ，类似于上面的推导可得 n n
n n e e e e ??
?
??-=-=-51,5101
每递推一步误差缩小到原值的
5
1
，所以递推公式(**)是数值稳定的。 由于]1,0[∈x 时，
5
56n
n n x x x x <+< 所以有估计式
)
1(51
)1(61+<<+n I n n
于是
21
51
216120?<
0087301587
.02
1051126120≈+
≈I
可得另一算法：??
?
??=-=≈-)1,,19,20()1
(510087301587.0120 n I n I I n n 这是数值稳定性好的算法。
3. 用秦九韶法计算123456)(23455+++++=x x x x x x P 在5=x 处的值，并验证之。 答案：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. 已知函数x y 2=的函数表
（1）试以02=x ，13=x 为节点，建立线性插值多项式)(1x L ，用以计算3
.02的近似
值，并估计截断误差。
（2）以11-=x ，02=x ，13=x 为节点，建立抛物线插值多项式)(2x L ，用以计算3
.02的近似值，并估计截断误差。
解（1）线性插值多项式为
120
10
1101)(1+=?--+?--=
x x x x L 则
3.1)3.0(213.0=≈L
其截断误差由(2-13)得
))((2
)(322
1x x x x M x R --≤
其中)("max 3
22x f M x x x ≤≤=，因为x x f 2)2(ln )("2=，可得
9069.0)2(ln 2)("max 21
0==≤≤x f x
于是
09522.0)13.0(3.02
9069
.0)3.0(1=-≤
R （2）抛物插值多项式为
175.025.02)
11)(01()
1)(0(1)10)(10()1)(1(5.0)11)(01()1)(0()(22++=?+-+-+?-+-++?------=
x x x x x x x x x L
则
248.1)3.0(223.0=≈L
其截断误差为
))()((6
)(3213
2x x x x x x M x R ---≤
其中6660.0)2(ln 2)('''max 3
1
13===≤≤-x f M x ，
故有
03030.0)13.0)(03.0)(13.0(6
6660
.0)3.0(2=--+≤
R
5. 已知函数表 求出Lagrange 插值多项式，并计算2.1=x 处的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) 插值多项式的简式为
3
20833333.00833333.083333.46667.20x x x ++-
插值为 15.1307
6. 证明∑
5
=20=)()(i i i x l x x －，其中)(x l i )5,,1,0=( i 为插值基函数。 证明
∑∑
5
=225
=2
)()+2(=)()(i i i i i i i x l x xx x x l x x －－
0+2=)(+)(2)(225
=5
=5
=2
2
∑
∑
∑＝－－＝x xx x x l x x l x x x l x i i i i i i i i
7. 设],[)(4b a C x f ∈，已知插值节点 b x x a n <<<<0 且3,2,1=,1i h x x i i ＝－－，
证明：（1）)(''max 8
)(≤
)()(max ≤≤2
1≤≤x f a b x L x f b x a b x a －－ （2）)('''max 27
3≤
)()(max 202
0≤≤3
2≤≤x f h x L x f x x x x x x － （3）)(max 241≤)()(max )4(≤≤4
3≤≤2
030x f h x L x f x x x x x x －
证明 （1）由插值余项得
))((!
2)
(''=
)()(=)(11b x a x ξf x L x f x R －－－ b ξa << 又 4
)(=))((=))((max 2
2
+=≤≤a b b x a x b x a x b
a x b
x a －－－－－
故 )(''m a x 8)(≤
)(''m a x ))((max 21
≤)()(max ≤≤2≤≤≤≤1≤≤x f a b x f b x a x x L x f b x a b x a b x a b x a －－－－
当x x f =)(时，等式成立。
（2）由插值余项得
))()((!
3)
('''=
)()(=)(21022x x x x x x ξf x L x f x R －－－－ 20< )('''max )2)((max 61 ≤ )()(max 202 02 0≤≤≤≤2≤≤x f h y h y y x L x f x x x x x x x x x －－－ 而 )2)((=)(h y h y y y g －－在h y 2≤≤0上的最大值为 3 9 32h 从而有 )('''m a x )2)((max 61 ≤ )()(max 202 02 0≤≤≤≤2≤≤x f h y h y y x L x f x x x x x x x x x －－－ ＝ 3 9 32h )('''max 2 0≤≤x f x x x （3）由插值余项得 ))()()((! 4) (= )()(=)(3210) 4(33x x x x x x x x ξf x L x f x R －－－－－ 30< )(max )3)(2)((max !41 ≤ )()(max )4(≤≤≤≤3≤≤303 03 0x f h y h y h y y x L x f x x x x x x x x x －－－－ =)(max )2+3)(3(max !41 )4(≤≤222≤≤303 0x f h hy y hy y x x x x x x －－ = )(max )2+(max ! 41 )4(≤≤20≤≤h 4 9 3 02x f h u u x x x u － 其中hy y u 32-= =)(max !41)4(≤≤4 3 0x f h x x x 8. 已知19个点的坐标如下表，根据表中数据求一插值多项式，并由此计算75.11=x 点的函数值的近似值。 数据表： 10 15 20 25 30 35 40 2.30259 2.70805 2.99573 3.21888 3.4012 3.55535 3.6888 45 50 55 60 65 70 75 3.80666 3.91202 4.00733 4.09434 4.17439 4.2485 4.31749 80 85 90 95 100 4.38203 4.44265 4.49981 4.55388 4.60517 解 已知数据表： 10 15 20 25 30 35 40 2.30259 2.70805 2.99573 3.21888 3.4012 3.55535 3.6888 45 50 55 60 65 70 75 3.80666 3.91202 4.00733 4.09434 4.17439 4.2485 4.31749 80 85 90 95 100 4.38203 4.44265 4.49981 4.55388 4.60517 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. 所求插值445.2)75.11(≈f 。 9. 已知函数x y cos =的函数表如下，试用4次等距节点插值公式计算)048 .0(f 的近似值，并估计误差。 解 先根据表2-3构造差分表。 用Newton 向前插值公式计算)048.0(f 的近似值，取048.0=x ，1.0=h ， 48.00 =-= h x t ，用差分表的上半部分差分，得 048 .0cos 99885.000012 .0! 4) 348.0)(248.0)(148.0(48.000013 .0!3)248.0)(148.0(48.0) 00993.0(2 ) 148.0(48.0)00500.0(48.00000.1)048.0(4≈=?---?+ ?--?+-?-?+-?+=N 误差估计由式(2-38)可得 755 4105487.1)4)(3)(2)(1(! 5)048.0(-?≤----≤ h t t t t t M R 其中565.06.0sin 5≤=M 。 用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 化简得： 64320000.000117.000354.000526.000279.01t t t t t ++--+ 1.00012 0.99942 即 ≈)048 .0(f 1.00012。 综合实习题（阅读） （1991年美国大学生数学建模竞赛A 题） 某些州的用水管理机构需估计公众的用水速度（单位：G/h ）和每天总用水量。许多供 水单位由于没有测量流入或流出水箱流量的设备，而只能测量水箱中的水位（误差不超过5％）。当水箱水位低于最低水位L 时，水泵开始工作将水灌入水箱，直至水位达到最高水位H 为止。但是依然无法测量水泵灌水流量，因此，在水泵工作时无法立即将水箱中的水位和水量联系起来。水泵一天灌水1~2次，每次约两小时。试估计在任一时刻（包括水泵灌水期间）t 流出水箱的流量)(t f ，并估计一天的总用水量。 【注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 －水的容积 i V －时刻i t （h ）水的容积 )(t f －流出水箱的水的流速（G/h ） p －水泵的灌水速度（G/h ） 为讨论方便，将表1中数据的单位变换为： 时间单位：h 水量： h r V 2π=，单位：G 310 变换结果见表2 表2 为求函数)(t f ，先求相邻时间区间的中点27,..2,121=+=+i t t t i i i 及在时间区间内水箱中流出的水的平均速度 27,..2,111=--=++i t t v v v i i i 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 表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 2t + 217.299 3t - 54.55444t + 6.389795t - 0.379256 6t + 0.0109854 7t - 0.000123258 8t 由图可见，将460556.01=t 小时和460556.242=t 小时代入)(t f 得 13339)(1=t f ，3.10933)(2=t f 3.2395)()(12=-t f t f 一天内用水总量近似等于 ? =460556 .24460556 .0)(dt t f Q ＝275006 Integrate[ft,{t,0.460556,24.460556}]//N 若按常规每千人的用水量为100000G/h 计算，因此估计该地区大约有2780人。 目录 第一章非线性方程求根 (3) 1.1迭代法 (3) 1.2牛顿法 (4) 1.3弦截法 (5) 1.4二分法 (6) 第二章插值 (7) 2.1线性插值 (7) 2.2二次插值 (8) 2.3拉格朗日插值 (9) 2.4分段线性插值 (10) 2.5分段二次插值 (11) 第三章数值积分 (13) 3.1复化矩形积分法 (13) 3.2复化梯形积分法 (14) 3.3辛普森积分法 (15) 3.4变步长梯形积分法 (16) 第四章线性方程组数值法 (17) 4.1约当消去法 (17) 4.2高斯消去法 (18) 4.3三角分解法 (20) 4.4雅可比迭代法 (21) 4.5高斯—赛德尔迭代法 (23) 第五章常积分方程数值法 (25) 5.1显示欧拉公式法 (25) 5.2欧拉公式预测校正法 (26) 5.3改进欧拉公式法 (27) 5.4四阶龙格—库塔法 (28) 数值计算方法 第一章非线性方程求根 1.1迭代法 程序代码： Private Sub Command1_Click() x0 = Val(InputBox("请输入初始值x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = (Exp(2 * x0) - x0) / 5 If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求f(x)=e2x-6x=0在x=0.5附近的根（ep=10-10） 1.2牛顿法 程序代码： Private Sub Command1_Click() b = Val(InputBox("请输入被开方数x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = x0 - (x0 ^ 2 - b) / (2 * b) If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求56的值。（ep=10-10） 大作业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用 程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) 《数值计算方法》复习试题 一、填空题： 1、????? ?????----=410141014A ，则A 的LU 分解为 A ??? ?????????=? ?????????? ?。 答案： ?? ????????--??????????--=1556141501 4115401411A 2、已知3.1)3(,2.1)2(,0.1)1(===f f f ，则用辛普生（辛卜生）公式计算求得 ?≈3 1 _________ )(dx x f ,用三点式求得≈')1(f 。 答案：， 3、1)3(,2)2(,1)1(==-=f f f ，则过这三点的二次插值多项式中2 x 的系数为 ， 拉格朗日插值多项式为 。 答案：-1， )2)(1(21 )3)(1(2)3)(2(21)(2--------= x x x x x x x L 4、近似值*0.231x =关于真值229.0=x 有( 2 )位有效数字； 5、设)(x f 可微,求方程)(x f x =的牛顿迭代格式是( )； （ 答案 )(1)(1n n n n n x f x f x x x '--- =+ 6、对1)(3 ++=x x x f ,差商=]3,2,1,0[f ( 1 ),=]4,3,2,1,0[f ( 0 )； 7、计算方法主要研究( 截断 )误差和( 舍入 )误差； 8、用二分法求非线性方程 f (x )=0在区间(a ,b )内的根时，二分n 次后的误差限为 ( 1 2+-n a b )； 9、求解一阶常微分方程初值问题y '= f (x ,y )，y (x 0)=y 0的改进的欧拉公式为 ( )] ,(),([2111+++++=n n n n n n y x f y x f h y y )； 10、已知f (1)＝2，f (2)＝3，f (4)＝，则二次Newton 插值多项式中x 2系数为( )； 11、 两点式高斯型求积公式?1 d )(x x f ≈( ?++-≈1 )] 321 3()3213([21d )(f f x x f )，代数精 度为( 5 )； 12、 解线性方程组A x =b 的高斯顺序消元法满足的充要条件为(A 的各阶顺序主子式均 不为零)。 13、 为了使计算 32)1(6 )1(41310-- -+-+ =x x x y 的乘除法次数尽量地少，应将该表 达式改写为 11 ,))64(3(10-= -++=x t t t t y ，为了减少舍入误差，应将表达式 19992001-改写为 199920012 + 。 14、 用二分法求方程01)(3 =-+=x x x f 在区间[0,1]内的根,进行一步后根的所在区间 为 ，1 ,进行两步后根的所在区间为 ， 。 15、 、 16、 计算积分?1 5 .0d x x ,取4位有效数字。用梯形公式计算求得的近似值为 ，用辛卜 生公式计算求得的近似值为 ，梯形公式的代数精度为 1 ，辛卜生公式的代数精度为 3 。 17、 求解方程组?? ?=+=+042.01532121x x x x 的高斯—塞德尔迭代格式为 ?????-=-=+++20/3/)51()1(1)1(2)(2)1(1 k k k k x x x x ，该迭 代格式的迭代矩阵的谱半径)(M ρ= 121 。 18、 设46)2(,16)1(,0)0(===f f f ,则=)(1x l )2()(1--=x x x l ，)(x f 的二次牛顿 插值多项式为 )1(716)(2-+=x x x x N 。 19、 求积公式 ?∑=≈b a k n k k x f A x x f )(d )(0 的代数精度以( 高斯型 )求积公式为最高，具 有( 12+n )次代数精度。 大作业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) while flag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1; while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +?? =-= ?-???解： Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; 数值计算方法作业 姓名：李琦 学号：062410124 求 013=--x x 在x=1.5附近的一个根。 一．牛顿下山法： #include 二．加权法 #include 三．单点弦法： #include 参考答案 第一章 1 *1x =1.7； * 2x =1.73； *3x =1.732 。 2． 3． （1） ≤++)(* 3*2*1x x x e r 0.00050； （注意：应该用相对误差的定义去求） （2） ≤)(*3*2*1x x x e r 0.50517； （3） ≤)/(*4*2x x e r 0.50002。 4．设6有n 位有效数字，由6≈2.4494……，知6的第一位有效数字1a =2。 令3)1()1(1* 102 1 102211021)(-----?≤??=?= n n r a x ε 可求得满足上述不等式的最小正整数n =4，即至少取四位有效数字，故满足精度要求可取6≈2.449。 5. 答：（1）*x (0>x )的相对误差约是* x 的相对误差的1/2倍； （2）n x )(* 的相对误差约是* x 的相对误差的n 倍。 6． 根据******************** sin 21)(cos 21sin 21)(sin 21sin 21)(sin 21)(c b a c e c b a c b a b e c a c b a a e c b S e r ++≤ =* *****) ()()(tgc c e b b e a a e ++ 注意当20* π < 7．设20= y ，41.1*0 =y ，δ=?≤--2* 00102 1y y 由 δ1* 001*111010--≤-=-y y y y ， δ2*111*221010--≤-=-y y y y M δ10*991*10101010--≤-=-y y y y 即当0y 有初始误差δ时，10y 的绝对误差的绝对值将减小10 10-倍。而110 10 <<-δ，故计算过程稳定。 8. 变形后的表达式为： （1）)1ln(2--x x =)1ln(2-+-x x （2）arctgx x arctg -+)1(=) 1(11 ++x x arctg （3） 1ln )1ln()1(ln 1 --++=? +N N N N dx x N N =ΛΛ+-+- +3 2413121)1ln(N N N N 1ln )11ln()1(-++ +=N N N N =1)1ln()1 1ln(-+++N N N （4）x x sin cos 1-=x x cos 1sin +=2x tg 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2; 《数值计算方法》复习试题 一、填空题: 1、????? ?????----=410141014A ,则A 的LU 分解为 A ??? ?????????=? ?????????? ?。 答案: ?? ????????--??????????--=1556141501 4115401411A 3、1)3(,2)2(,1)1(==-=f f f ,则过这三点的二次插值多项式中2 x 的系数为 ,拉 格朗日插值多项式为 。 答案:-1, )2)(1(21 )3)(1(2)3)(2(21)(2--------= x x x x x x x L 4、近似值*0.231x =关于真值229.0=x 有( 2 )位有效数字; 5、设)(x f 可微,求方程)(x f x =的牛顿迭代格式就是( ); 答案 )(1)(1n n n n n x f x f x x x '--- =+ 6、对1)(3 ++=x x x f ,差商=]3,2,1,0[f ( 1 ),=]4,3,2,1,0[f ( 0 ); 7、计算方法主要研究( 截断 )误差与( 舍入 )误差; 8、用二分法求非线性方程 f (x )=0在区间(a ,b )内的根时,二分n 次后的误差限为 ( 1 2+-n a b ); 10、已知f (1)＝2,f (2)＝3,f (4)＝5、9,则二次Newton 插值多项式中x 2系数为( 0、15 ); 11、 解线性方程组A x =b 的高斯顺序消元法满足的充要条件为(A 的各阶顺序主子式均 不为零)。 12、 为了使计算 32)1(6 )1(41310-- -+-+ =x x x y 的乘除法次数尽量地少,应将该表 达式改写为 11 ,))64(3(10-= -++=x t t t t y ,为了减少舍入误差,应将表达式 2020年奥鹏吉大网络教育《计算方法》大作业解答 (说明：前面是题目，后面几页是答案完整解答部分，注意的顺序。) 一、解线性方程 用矩阵的LU分解算法求解线性方程组 用矩阵的Doolittle分解算法求解线性方程组 用矩阵的Doolittle分解算法求解线性方程组 用高斯消去法求解线性方程组 用高斯消去法求解线性方程组 用主元素消元法求解线性方程组 用高斯消去法求解线性方程组 利用Doolittle分解法解方程组Ax=b，即解方程组 1、用矩阵的LU分解算法求解线性方程组 X1+2X2+3X3 = 0 2X1+2X2+8X3 = -4 -3X1-10X2-2X3 = -11 2、用矩阵的Doolittle分解算法求解线性方程组 X1+2X2+3X3 = 1 2X1– X2+9X3 = 0 -3X1+ 4X2+9X3 = 1 3、用矩阵的Doolittle分解算法求解线性方程组 2X1+X2+X3 = 4 6X1+4X2+5X3 =15 4X1+3X2+6X3 = 13 4、用高斯消去法求解线性方程组 2X 1- X 2+3X 3 = 2 4X 1+2X 2+5X 3 = 4 -3X 1+4X 2-3X 3 = -3 5、用无回代过程消元法求解线性方程组 2X 1- X 2+3X 3 = 2 4X 1+2X 2+5X 3 = 4 -3X 1+4X 2-3X 3 = -3 6、用主元素消元法求解线性方程组 2X 1- X 2+3X 3 = 2 4X 1+2X 2+5X 3 = 4 -3X 1+4X 2-3X 3 = -3 7、用高斯消去法求解线性方程组 123123123234 4272266 x x x x x x x x x -+=++=-++= 8、利用Doolittle 分解法解方程组Ax=b ，即解方程组 12341231521917334319174262113x x x x -? ????? ???? ??-??????=? ? ????--?????? --???? ?? 第二章 1. 题目：运用MATLAB编程实现牛顿迭代 2. 实验操作 1、打开MATLAB程序软件。 2、在MATLAB中编辑如下的M程序。 function [p1,err,k,y]=newton(f,df,p0,delta,max) %f 是要求根的方程(f(x)=0); %df 是f(x)的导数； %p0是所给初值,位于x*附近； %delta是给定允许误差； %max是迭代的最大次数； %p1是newton法求得的方程的近似解； %err是p0的误差估计； %k是迭代次数； p0 for k=1:max p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1; k p1 err y=feval('f',p1) if (err p0 = 1.2000 k =1 p1=1.1030 err=0.0970 y=0.0329 k= 2 p1=1.0524 err=0.0507 y=0.0084 k =3 p1=1.0264 err=0.0260 y=0.0021 k =4 p1=1.0133 err=0.0131 y=5.2963e-004 k =5 p1=1.0066 err=0.0066 y=1.3270e-004 k =6 p1=1.0033 err=0.0033 y=3.3211e-005 k =7 p1=1.0017 err=0.0017 y=8.3074e-006 k =8 p1=1.0008 err=8.3157e-004 y = 2.0774e-006 k =9 p1=1.0004 err=4.1596e-004 y =5.1943e-007 k=10 p1=1.0002 err=2.0802e-004 y= 1.2987e-007 k=11 p1=1.0001 err=1.0402e-004 y =3.2468e-008 k=12 p1=1.0001 err=5.2014e-005 y=8.1170e-009 k=13 p1=1.0000 err=2.6008e-005 y= 2.0293e-009 k=14 p1=1.0000 err=1.3004e-005 y=5.0732e-010 k=15 p1 =1.0000 err=6.5020e-006 y=1.2683e-010 k=16 p1 =1.0000 err=3.2510e-006 y=3.1708e-011 k=17 p1 =1.0000 err=1.6255e-006 y =7.9272e-012 k=18 p1 =1.0000 err =8.1279e-007 y= 1.9820e-012 ans = 1.0000 结果说明：经过18次迭代得到精确解为1，误差为8.1279e-007。 数值分析报大作业 班级：铁道2班 专业：道路与铁道工程 姓名：蔡敦锦 学号：13011260 一、序言 该数值分析大作业是通过C语言程序编程在Microsoft Visual C++ 6.0编程软件上运行实现的。本来是打算用Matlab软间来计算非线性方程的根的。学习Matlab也差不多有一个多月了，感觉自己编程做题应该没什么问题了；但是当自己真心的去编程、运行时才发现有很多错误，花了一天时间修改、调试程序都没能得到自己满意的结果。所以，我选择了自己比较熟悉的C程序语言来编程解决非线性的求值问题，由于本作业是为了比较几种方法求值问题的收敛速度和精度的差异，选择了一个相对常见的非线性函数来反映其差异，程序运行所得结果我个人比较满意。编写C语言，感觉比较上手，程序出现问题也能比较熟练的解决。最终就决定上交一份C程序语言编程的求值程序了！ 二、选题 本作业的目的是为了加深对非线性方程求根方法的二分法、简单迭代法、、牛顿迭代法弦截法等的构造过程的理解；能将各种方法的算法描述正确并且能够改编为程序并在计算机上实现程序的正确合理的运行，能得到自己满意的结果，并且能调试修改程序中可能出现的问题和程序功能的增减修改。本次程序是为了比较各种方法在求解同一非线性方程根时，在收敛情况上的差异。 为了达到上面的条件我选择自己比较熟悉的语言—C语言来编程，所选题目为计算方程f(x)=x3-2x-5=0在区间[2,3]内其最后两近似值的差的绝对值小于等于5 ?的根的几种方法的比较。 110- 本文将二分法、牛顿法、简单迭代法、弦截法及加速收敛法这五种方法在同一个程序中以函数调用的方式来实现，比较简洁明了，所得结果能很好的比较，便于分析；发现问题和得出结论。 《计算方法》期中复习试题 一、填空题： 1、已知3.1)3(,2.1)2(,0.1)1(===f f f ，则用辛普生（辛卜生）公式计算求得 ?≈3 1 _________ )(dx x f ,用三点式求得≈')1(f 。 答案：， 2、1)3(,2)2(,1)1(==-=f f f ，则过这三点的二次插值多项式中2 x 的系数为 ， 拉格朗日插值多项式为 。 答案：-1， )2)(1(21 )3)(1(2)3)(2(21)(2--------= x x x x x x x L 3、近似值*0.231x =关于真值229.0=x 有( 2 )位有效数字； 4、设)(x f 可微,求方程)(x f x =的牛顿迭代格式是( )； 答案 )(1)(1n n n n n x f x f x x x '--- =+ 5、对1)(3 ++=x x x f ,差商=]3,2,1,0[f ( 1 ),=]4,3,2,1,0[f ( 0 )； 6、计算方法主要研究( 截断 )误差和( 舍入 )误差； 7、用二分法求非线性方程f (x )=0在区间(a ,b )内的根时，二分n 次后的误差限为 ( 1 2+-n a b )； 8、已知f (1)＝2，f (2)＝3，f (4)＝，则二次Newton 插值多项式中x 2系数为( )； 11、 两点式高斯型求积公式?1 d )(x x f ≈( ?++-≈1 )] 321 3()3213([21d )(f f x x f )，代数精 度为( 5 )； 12、 为了使计算 32)1(6 )1(41310-- -+-+ =x x x y 的乘除法次数尽量地少，应将该表 达式改写为 11 ,))64(3(10-= -++=x t t t t y ，为了减少舍入误差，应将表达式 19992001-改写为 199920012 + 。 计算方法上机报告 姓名： 学号： 班级： 目录 题目一------------------------------------------------------------------------------------------ - 4 - 1.1题目内容 ---------------------------------------------------------------------------- - 4 - 1.2算法思想 ---------------------------------------------------------------------------- - 4 - 1.3Matlab源程序----------------------------------------------------------------------- - 5 - 1.4计算结果及总结 ------------------------------------------------------------------- - 5 - 题目二------------------------------------------------------------------------------------------ - 7 - 2.1题目内容 ---------------------------------------------------------------------------- - 7 - 2.2算法思想 ---------------------------------------------------------------------------- - 7 - 2.3 Matlab源程序---------------------------------------------------------------------- - 8 - 2.4计算结果及总结 ------------------------------------------------------------------- - 9 - 题目三----------------------------------------------------------------------------------------- - 11 - 3.1题目内容 --------------------------------------------------------------------------- - 11 - 3.2算法思想 --------------------------------------------------------------------------- - 11 - 3.3Matlab源程序---------------------------------------------------------------------- - 13 - 3.4计算结果及总结 ------------------------------------------------------------------ - 14 - 题目四----------------------------------------------------------------------------------------- - 15 - 4.1题目内容 --------------------------------------------------------------------------- - 15 - 4.2算法思想 --------------------------------------------------------------------------- - 15 - 4.3Matlab源程序---------------------------------------------------------------------- - 15 - 4.4计算结果及总结 ------------------------------------------------------------------ - 16 - 题目五----------------------------------------------------------------------------------------- - 18 - 第四章 问题一 一、问题综述 在离地球表面高度为y处的重力加速度如下： 计算高度y=55000m处的重力加速度值。 二、问题分析 以高度y作为自变量，重力加速度的值为因变量。得到以下信息： f(0)=9.8100; f(30000)=9.7487; f(60000)=9.6879; f(90000)=9.6278; f(120000)=9.5682; 本题要求的就是f(55000)的值。 以下将采用课堂中学到的Lagrange插值多项式法、Newton插值多项式法、分段低次插值法和样条插值法求解该问题。 三、问题解决 1. lagrange插值多项式法 对某个多项式函数，已知有给定的k+ 1个取值点： 其中对应着自变量的位置，而对应着函数在这个位置的取值。 假设任意两个不同的x j都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为： 其中每个为拉格朗日基本多项式（或称插值基函数），其表达式为： 拉格朗日基本多项式的特点是在上取值为1，在其它的点上取值为0。 源程序lagrange.m function [c,f]=lagrange(x,y,a) % 输入：x是自变量的矩阵；y是因变量的矩阵；a是要计算的值的自变量； % 输出：c是插值多项式系数矩阵；f是所求自变量对应的因变量； m=length(x); l=zeros(m,m); % l是权矩阵 f=0; for i=1:m v=1; for j=1:m if i~=j v=conv(v,poly(x(j)))/(x(i)-x(j)); % v是l_i(x)的系数矩阵 end end l(i,:)=v; % l矩阵的每一行都是x从高次到低次的系数矩阵 end c=vpa(y*l,10); % 对应阶次的系数相加，乘以y，显示10位有效数字 for k=1:m f=f+c(k)*a^(m-k); end 输入矩阵 x=[0 30000 60000 90000 120000] y=[9.81 9.7487 9.6879 9.6278 9.5682] a=55000 再运行源函数，可得： c = [ -2.057613169e-23, 4.938271605e-18, -3.703703702e-14, -0.000002046111111, 9.81] f = 9.6979851723251649906109417384537 计算方法实验报告 班级： 学号： 姓名： 成绩： 1 舍入误差及稳定性 一、实验目的 （1）通过上机编程，复习巩固以前所学程序设计语言及上机操作指令； （2）通过上机计算，了解舍入误差所引起的数值不稳定性 二、实验内容 1、用两种不同的顺序计算10000 21n n -=∑，分析其误差的变化 2、已知连分数() 1 01223//(.../)n n a f b b a b a a b =+ +++，利用下面的算法计算f ： 1 1 ,i n n i i i a d b d b d ++==+ (1,2,...,0 i n n =-- 0f d = 写一程序，读入011,,,...,,,...,,n n n b b b a a 计算并打印f 3、给出一个有效的算法和一个无效的算法计算积分 1 041 n n x y dx x =+? (0,1,...,1 n = 4、设2 2 11N N j S j == -∑ ，已知其精确值为1311221N N ?? -- ?+?? （1）编制按从大到小的顺序计算N S 的程序 （2）编制按从小到大的顺序计算N S 的程序 （3）按两种顺序分别计算10001000030000,,,S S S 并指出有效位数 三、实验步骤、程序设计、实验结果及分析 1、用两种不同的顺序计算10000 2 1n n -=∑，分析其误差的变化 （1）实验步骤： 分别从1~10000和从10000~1两种顺序进行计算，应包含的头文件有stdio.h 和math.h （2）程序设计： a.顺序计算 #include 数值计算大作业 题目一、非线性方程求根 1.题目 假设人口随时间和当时人口数目成比例连续增长，在此假设下人口在短期内的增长建立数学模型。 （1）如果令()N t 表示在t 时刻的人口数目，β 表示固定的人口出生率，则人口数目满足微分方程() ()dN t N t dt β=，此方程的解为0()=t N t N e β； （2）如果允许移民移入且速率为恒定的v ，则微分方程变成() ()dN t N t v dt β=+， 此方程的解为 0()=+ (1) t t v N t N e e βββ -； 假设某地区初始有1000000人，在第一年有435000人移入，又假设在第一年年底该地区人口数量1564000人，试通过下面的方程确定人口出生率β，精确到 410-；且通过这个数值来预测第二年年末的人口数，假设移民速度v 保持不变。 435000 1564000=1000000(1) e e βββ + - 2.数学原理 采用牛顿迭代法，牛顿迭代法的数学原理是，对于方程0)(=x f ，如果) (x f 是线性函数，则它的求根是很容易的，牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程0)(=x f 逐步归结为某种线性方程来求解。 设已知方程0)(=x f 有近似根k x (假定0)(≠'x f )，将函数)(x f 在点k x 进行泰勒展开，有 . ))(()()(???+-'+≈k k k x x x f x f x f 于是方程0)(=x f 可近似地表示为 ))(()(=-'+k k x x x f x f 这是个线性方程，记其根为1k x +，则1k x +的计算公式为 数值分析复习试题 第一章 绪论 一. 填空题 1.* x 为精确值 x 的近似值；() **x f y =为一元函数 ()x f y =1的近似值； ()**,*y x f y =为二元函数()y x f y ,2=的近似值，请写出下面的公式：**e x x =-： *** r x x e x -= ()()()*'1**y f x x εε≈? ()() () ()'***1**r r x f x y x f x εε≈ ? ()()()() ()* *,**,*2**f x y f x y y x y x y εεε??≈?+??? ()()()()() ** * *,***,**222r f x y e x f x y e y y x y y y ε??≈ ?+??? 2、 计算方法实际计算时，对数据只能取有限位表示，这时所产生的误差叫 舍入误 差 。 3、 分别用2.718281，2.718282作数e 的近似值，则其有效数字分别有 6 位和 7 位；又取 1.73≈-21 1.73 10 2 ≤?。 4、 设121.216, 3.654x x ==均具有3位有效数字，则12x x 的相对误差限为 0.0055 。 5、 设121.216, 3.654x x ==均具有3位有效数字，则12x x +的误差限为 0.01 。 6、 已知近似值 2.4560A x =是由真值T x 经四舍五入得 到,则相对误差限为 0.0000204 . 7、 递推公式,??? ? ?0n n-1y =y =10y -1,n =1,2, 如果取0 1.41y ≈作计算,则计算到10y 时,误 差为 81 10 2 ?;这个计算公式数值稳定不稳定 不稳定 . 8、 精确值 14159265.3* =π，则近似值141.3*1=π和1415.3*2=π分别有 3 计算方法大作业 题目：非线性方程求根的新方法 班级：xxx 学号：xxx 姓名：xxx 非线性方程求根的新方法 一、问题引入 在计算和实际问题中经常遇到如下非线性问题的求解： F(x)=0 (1) 我们经常采用的方法是经典迭代法： 经典迭代方法 不动点迭代方法是一种应用广泛的方法，其加速方法较多，如Stiffensen加速方法的局部收敛阶(以下简称为收敛阶)为2阶；牛顿迭代方法的收敛阶亦为2阶，且与其相联系的一些方法如简化牛顿法、牛顿下山法、弦截法的收敛阶阶数介于1和2之间；而密勒法的收敛阶与牛顿法接近，但计算量较大且涉及零点的选择问题，同时收敛阶也不够理想。 因此本文介绍一种新的迭代方法 从代数角度看，牛顿法和密勒法分别是将f(x)在xk附近近似为一线性函数和二次抛物插值函数，一种很自然的想法就是能否利用Taylor展开，将f(x)在xk附近近似为其他的二次函数?答案是肯定的.其中的一种方法是将f(x)在Xk处展开3项，此时收敛阶应高于牛顿法，这正是本文的出发点． 二、算法推导 设函数f(x)在xk附近具有二阶连续导数，则可将f(x)在xk处进行二阶Taylor展开，方程(1) 可近似为如下二次方程： f(xk)+f’(xk)(x-xk)+2^(-1)f’’(xk)(x-xk)^2=0，(2) 即 2^(-1)f’’(xk)x^2+(f’(xk)-xkf’’(xk))x+2^(-1)f’’(xk)xk^2-xkf’(xk)+f(xk)=0(3) 利用求根公式可得 X=xk-(f’’(xk))^(-1)(f’(xk))-sqrt((f’(xk)^2±2f’’(xk)f(xk)))(4) 其中±符号的选取视具体问题而定，从而可构造迭代公式 X k+1=xk-(f’’(xk))^(-1)(f’(xk))-sqrt((f’(xk)^2±2f’’(xk)f(xk)))(5) 确定了根号前正负号的迭代公式(5)，可称为基于牛顿法和Taylor展开的方法，简记为BNT 方法． 为描述方便起见，以下将f(xk)，f’(xk)，f’’(xk)分别记为f,f’,f’’．首先，二次方程(3)对应于一条抛物曲线，其开口方向由f’’(xk)，x∈U(xk)的符号确定，其中U(xk)为xk的某邻域，其顶点为 P(xk-(f’’)^(-1)f’，fk-(2f’’)^(-1)(f’)^2)．为使(5)式唯一确定x k+1，须讨论根式前正负号的取舍问题．下面从该方法的几何意义分析(5)式中正负号的取舍． 1)当f(xk)=o时，z。即为所求的根． 2)当f(xk)>O时，根据y=f(x)的如下4种不同情形(见图1)确定(5)式中根号前的符号． (a)当f’’(xk)数值计算方法大作业
数值分析大作业-三、四、五、六、七
《数值计算方法》试题集及答案
数值分析大作业三 四 五 六 七
数值计算方法作业
西工大计算方法作业答案
(完整版)数值计算方法上机实习题答案
《数值计算方法》试题集及答案
2020年奥鹏吉大网络教育《计算方法》大作业解答
数值分析作业
数值分析大作业
数值计算方法试题集和答案
西安交通大学计算方法B大作业
数值计算方法第4次作业
计算方法上机实习题大作业(实验报告).
数值计算大作业
数值分析计算方法试题集及答案
计算方法大作业非线性方程求根的新方法
相关主题
文本预览