研究生“数值分析”课后题(上机编程部分)答案
数值分析课后习题答案
0 1
0 10 1 1 0 0 0 1
0 0 12 1 1 2 0 0 0
1 2
0 0 0 1 1 0
1 2
1 2
1 2
1
0 0 0 1 0
1 2
1 2
0
1 2
1 2
0
0
0
341 1 1
2-5.对矩阵A进行LDLT分解和GGT分解，并求解方程组
Ax=b，其中
16 4 8
1
A 4 5 4 , b 2
8 4 22
3
解
16 A 4
4 5
84
44 11
2-3(1).对矩阵A进行LU分解,并求解方程组Ax=b,其中
2 1 1 A1 3 2
4 ,b6
1 2 2
5
解
2 A 1
1 3
1 2
2 11
22
1
5 2
1
3 21来自，所以 A12
1
2 1 1
5 3
2-2(1).用列主元Gauss消元法解方程组
3 2 6x1 4 10 7 0x2 7 5 1 5x3 6
解
3 2 6 4 10 7 0 7 10 7 0 7
r1r2
消元
10 7 0 7 3 2 6 4 0 0.1 6 6.1
r=0.5101-n/3.162…<0.5101-n/3<0.01% 因此只需n=5.即取101/2=3.1623
研究生数值分析上机试题及解答
东华大学研究生数值分析试题（上机部分）A 卷2008年12月 时间：60分钟班级 学号 机号 姓名 得分 注意：要求写出M 函数（如果需要）、MATLAB 命令和计算结果。
1. 求下列方程组在0<α, β<1中的解⎩⎨⎧-=+=βαββααsin 2.0cos 7.0cos 2.0sin 7.0 命令fun=inline('[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))]','x'); [x,f,h]=fsolve(fun,[0.5 0.5]) 结果α=0.5265，β=0.50792命令>> fun=inline('c(1)+c(2)*x.^2','c','x'); >> x=[1.1 1.3 1.4 1.6 1.8]; >> y=[26 22 23 24 25];>> c=lsqcurvefit(fun,[0 0],x,y) 结果 c =23.7256 0.12873．求解下列微分方程组2(0)2013(0)1x x yx t y x yy '=-=⎧<<⎨'=+=⎩(结果只要求写出t =1时的解) 命令>> fun=inline('[y(1)-2*y(2);3*y(1)+y(2)]','t','y'); >> [t,y]=ode45(fun,[0 1], [2 1]) 结果x(1)=-5.6020, y(1)=2.15634．用定步长Gauss 积分法(课本123页)计算积分31e ln(1)x x dx -+⎰的近似值（等分数取4，每段取2个Gauss 点）。
命令fun=inline('exp(-x).*log(1+x)','x'); nagsint(fun,1,3,4,2) 结果 0.30865．矩阵改进平方根分解(课本25页)的计算公式为: d 1=a 11, 对i =2, 3, ⋯, n ,iki k ik ii i j ij ij j k jk ik ij ij l s a d i j d s l l s a s ∑∑-=-=-=-==-=1111,1,,2,1 ,/ ,试编写矩阵改进平方根分解的程序，并求矩阵1111551514A -⎛⎫ ⎪=-- ⎪ ⎪-⎝⎭的改进平方根分解。
数值分析课后部分习题答案
习题一（P.14）1. 下列各近似值均有4个有效数字,300.2,521.13,001428.0***===z y x ,试指出它们的绝对误差和相对误差限.解*20.001428=0.142810x -=⨯有4个有效数，即4n =，2m =-由有效数字与绝对误差的关系得绝对误差限为611101022m n --⨯=⨯， 由有效数字与相对误差的关系得相对误差限为(1)3111101022n a ---⨯=⨯； *213.521=0.1352110y =⨯有4个有效数，即4n =，2m =由有效数字与绝对误差的关系得绝对误差限为211101022m n --⨯=⨯， 由有效数字与相对误差的关系得相对误差限为(1)3111101022n a ---⨯=⨯； *12.300=0.230010z =⨯有4个有效数，即4n =，1m =由有效数字与绝对误差的关系得绝对误差限为311101022m n --⨯=⨯， 由有效数字与相对误差的关系得相对误差限为(1)3111101024n a ---⨯=⨯.2．下列各近似值的绝对误差限都是31021-⨯,试指出它们各有几位有效数字.***2.00021,0.032,0.00052x y z ===解*12.000210.20002110x ==⨯，即1m =由有效数字与绝对误差的关系得 311101022m n --⨯=⨯， 即3m n -=-，所以，2n =；*10.0320.3210y ==⨯，即1m =由有效数字与绝对误差的关系得 311101022m n --⨯=⨯， 即3m n -=-，所以，4n =；*30.000520.5210z -==⨯，即3m =-由有效数字与绝对误差的关系得 311101022m n --⨯=⨯， 即3m n -=-，所以，0n =.4.设有近似数35.2,84.1,41.2***===z y x 且都有3位有效数字，试计算***z y x S +=，问S 有几位有效数字.解 方法一因*1*1*12.41=0.24110, 1.840.18410, 2.350.23510x y z =⨯==⨯==⨯都有3位有效数字，即3n =，1m =，则211|(*)|101022m n e x --≤⨯=⨯，211|(*)|101022m n e y --≤⨯=⨯，211|(*)|101022m n e z --≤⨯=⨯，|(**)||*(*)*(*)|*|(*)|*|(*)|e y z z e y y e z z e y y e z ≈+≤+222112.3510 1.8410 2.0951022---≤⨯⨯+⨯⨯=⨯，221|(***)||(*)(**)|10 2.095102e x y z e x e y z --+≈+≤⨯+⨯1110.259510102--=⨯≤⨯， 又1***=2.41 1.84 2.350.673410x y z ++⨯=⨯，此时1m =，1m n -=-，从而得2n =.方法一因*1*1*12.41=0.24110, 1.840.18410, 2.350.23510x y z =⨯==⨯==⨯都有3位有效数字，即3n =，1m =，则211|(*)|101022m n e x --≤⨯=⨯，2110(*)2|(*)|=||* 2.41r e x e x x -⨯≤， 211|(*)|101022m n e y --≤⨯=⨯，2110(*)2|(*)|=||* 1.84r e y e y y -⨯≤，211|(*)|101022m n e z --≤⨯=⨯，2110(*)2|(*)|=||* 2.35r e z e z z -⨯≤|(**)||(*)(*)|r r r e y z e y e z ≈+，***|(***)||(*)(**)|******r r rx y z e x y z e x e y z x y z x y z +≈+++2.41 1.84 2.35|(*)||(*)+(*)|2.41 1.84 2.35 2.41 1.84 2.35r rr e x e y e z ⨯≤++⨯+⨯22211110 1.8410 2.35102222.41 1.84 2.35 2.41 1.84 2.35 2.41 1.84 2.35---⨯⨯⨯⨯⨯≤+++⨯+⨯+⨯20.385410-<⨯21102-<⨯，由有效数字与绝对误差的关系得2n =.5.序列{}n y 有递推公式),2,1(,1101 =-=-n y y n n若41.120≈=y （三位有效数字），问计算10y 的误差有多大，这个计算公式稳定吗？解 用0ε表示0y 的误差，由41.120≈=y ，得0=0.0042ε，由递推公式),2,1(,1101 =-=-n y y n n ，知计算10y 的误差为810=0.4210ε⨯，因为初始误差在计算的过程中被逐渐的放大，这个计算公式不稳定.习题2 ( P.84)3.证明()1nkk lx ==∑，对所有的x其中()k l x 为Lagrange 插值奇函数. 证明 令()1f x =，则()1i f x =， 从而 0()()()()nnn k k k k k L x l x f x l x ====∑∑，又(1)1()()()0(1)!n n n f R x x n ξω++==+，可得 ()()1n l x f x ==，从而()1nkk lx ==∑.4. 求出在=012x ，，和3处函数2()1f x x =+的插值多项式. 解 方法一 因为给出的节点个数为4，而2()1f x x =+从而余项(4)34()()()04!f R x x ξω==，于是233()()()()=+1L x f x R x f x x =-=（n 次插值多项式对次数小于或等于的多项式精确成立）.方法二 因为(0)1(1)2(2)5(3)10f f f f ====，，，， 而0(1)(2)(3)1()=-(1)(2)(3)(01)(02)(03)6x x x l x x x x ---=------，1(2)(3)1()=(2)(3)(10)(12)(13)2x x x l x x x x --=-----，2(1)(3)1()=-(1)(3)(20)(21)(23)2x x x l x x x x --=-----，3(1)(2)1()=(1)(2)(30)(31)(32)6x x x l x x x x --=-----，从而30123()()(0)()(1)()(2)()(3)L x l x f l x f l x f l x f =+++2=+1x .5. 设2()[,]f x C a b ∈且()()0f a f b ==，求证21max |()|()max |()|8a x ba xb f x b a f x ≤≤≤≤''≤-.证明 因()()0f a f b ==，则1()0L x =， 从而1()()()()()2!f f x R x x a x b ξ''==--，由极值知识得 21max |()|()max |()|8a x ba xb f x b a f x ≤≤≤≤''≤-6. 证明 (()())()()()(+)f x g x f x g x f x g x h ∆=⋅∆+∆⋅. 证明 由差分的定义(()())(+)()()()f xg x f xh g x h f x g x ∆=+-[(+)()()(+)][()()()()]f x h g x h f x g x h f x g x h f x g x =+-++-()()()(+)f x g x f x g x h =⋅∆+∆⋅或着 (()())(+)()()()f x g x f x h g x h f x g x ∆=+-[(+)()()()][()()()()]f x hg xh f x h g x f x h g x f x g x =+-+++- ()()()()f x h g x f x g x =+⋅∆+∆⋅7. 证明 n 阶差商有下列性质(a ) 如果()()F x cf x =，则0101[,,,][,,,]n n F x x x cf x x x =. (b ) 如果()()()F x f x g x =+，则010101[,,,][,,,][,,,]n n n F x x x f x x x g x x x =+.证明 由差商的定义 (a ) 如果()()F x cf x =，则12011010[,,,]-[,,,][,,,]n n n n F x x x F x x x F x x x x x -=-120110[,,,]-[,,,]n n n cf x x x cf x x x x x -=-120110[,,,]-[,,,]n n n f x x x f x x x c x x -=⋅-01[,,,]n cf x x x =.(b ) 如果()()()F x f x g x =+，则12011010[,,,]-[,,,][,,,]n n n n F x x x F x x x F x x x x x -=-12120110110[[,,,][,,,]]-[[,,,][,,,]]n n n n n f x x x g x x x f x x x g x x x x x --++=-12011120110,,,]-[,,,][,,,][,,,]+n n n n n n f x x x f x x x g x x x g x x x x x x x ---=--[ 0101[,,,][,,,]n n f x x x g x x x =+8. 设74()3431f x x x x =+++，求0172,2,,2]f [，0182,2,,2]f [.解 由P.35定理7的结论(2)，得7阶差商0172,2,,2]=3f [(()f x 的最高次方项的系数)，8阶差商0182,2,,2]=0f [(8阶以上的差商均等与0).9. 求一个次数不超过4次的多项式()P x ，使它满足：(0)(0)0P P '==，(1)(1)1P P '==，(2)1P =.解 方法一 先求满足插值条件(0)0P =，(1)=1P ，(2)1P =的二次插值多项式2()P x 213=22x -+(L-插值基函数或待定系数法)， 设()P x 22=()(1)(2)(1)(2)P x Ax x x Bx x x +--+--213=22x x -+2+(1)(2)(1)(2)Ax x x Bx x x --+-- 从而()P x '323=4B +(39)(641)(2)2x A B x A B x A -+-+-++，再由插值条件(0)0P '=，(1)1P '=，得3=,4A -1=,4B所以 ()P x 213=22x x -+231(1)(2)(1)(2)44x x x x x x ---+--， 即 ()P x 41=4x 332x -29+4x .方法二 设()P x 23401234=a a x a x a x a x ++++， 则 ()P x '231234=234a a x a x a x +++由插值条件(0)(0)0P P '==，(1)(1)1P P '==，(2)1P =，得010********0123400++++1+2+3+41+2+4+8+161a a a a a a a a a a a a a a a a =⎧⎪=⎪⎪=⎨⎪=⎪=⎪⎩ 解得 234931=,=-,=424a a a ， 从而()P x 41=4x 332x -29+4x . 方法三 利用埃尔米特插值基函数方法构造. 10. 下述函数()S x 在[1,3]上是3次样条函数吗？3232321,12()=92217,23x x x x S x x x x x ⎧-++≤≤⎨-+-+≤≤⎩ 解 因为22362,12()=31822,23x x x S x x x x ⎧-+≤≤'⎨-+-≤≤⎩， 66,12()=618,23x x S x x x -≤≤⎧''⎨-+≤≤⎩而12(2)=1=(2)S S ，12(2)=2=(2)S S ''，12(2)=6=(2)S S ''''， 又()S x 是三次函数，所以函数()S x 在[1,3]上是3次样条函数.补 设f (x )=x 4，试利用L-余项定理写出以-1,0,1,2为插值节点的三次插值多项式.解 因为 (4)34()()()(+1)(1)(2)4!f R x x x x x x ξω==--，从而3233()()()22L x f x R x x x x =-=+-习题3 ( P.159)1．设n k k x 0)}({=ϕ为],[b a 上具有权函数0)(≥x ω的正交多项式组且)(x k ϕ为首项系数为1的k 次的多项式，则n k k x 0)}({=ϕ于],[b a 线性无关.解 方法一 因为n k k x 0)}({=ϕ为],[b a 上具有权函数0)(≥x ω的正交多项式组,则其Gram 行列式不等于零，采用反证法：若{}n ϕϕϕ,,,10 于],[b a 线性相关，于是，存在不全为零,,,,10n c c c 使0011()()()0,[,]n n c x c x c x x a b ϕϕϕ+++=∈上式两边与i ϕ作内积得到0011(,)(,)(,)0(0,1,,)i i n i n c c c i n ϕϕϕϕϕϕ+++==，由于{}i c 不全为零，说明以上的齐次方程组有非零解),,,,(10n c c c 故系数矩阵的行列式为零，即{}0,,,10=n G ϕϕϕ 与假设矛盾.方法二 因为n k k x 0)}({=ϕ为],[b a 上具有权函数0)(≥x ω的正交多项式组,则其Gram 行列式不等于零，由( P.95)定理2得n k k x 0)}({=ϕ于],[b a 线性无关.2．选择α，使下述积分取得最小值1221()[],a x x dx α--⎰120()()x b e x dx α-⎰解1221()[]a x x dx αα-∂-∂⎰1221=[]x x dx αα-∂-∂⎰1221=2[]()x x x dx α--⋅-⎰5112=5x α-4=5α，令1221[]=0x x dx αα-∂-∂⎰，得=0α. 12()()x b e x dx αα∂-∂⎰120=()xe x dx αα∂-∂⎰1=2()()x e x x dx α-⋅-⎰2=23α- 令120()=0x e x dx αα∂-∂⎰，得=3α.3．设],3,1[,1)(∈=x xx f 试用},1{1x H 求)(x f 一次最佳平方逼近多项式.解 取权函数为()x x ω=(为了计算简便)，则32311(1,1)42x xdx ===⎰,33321126(1,)(,1)33x x x x dx ====⎰, 343311(,)204x x x x dx ===⎰,33111((),1)2f x xdx x x=⋅==⎰，3232111((),)42x f x x x dx x =⋅==⎰， 得法方程0126423264203a a ⎡⎤⎢⎥⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦⎢⎥⎣⎦，解得011211311a a ⎧=⎪⎪⎨⎪=-⎪⎩， 所以)(x f 的一次最佳平方逼近多项式1123()1111P x x =-. 8．什么常数C 能使得以下表达式最小？ ∑=-ni x i iCe x f 12))((解21(())i n x i i f x Ce C =∂-∂∑1=2(())()i i nx x i i f x Ce e =-⋅-∑， 令21(())=0i nx i i f x Ce C =∂-∂∑，得121()(),iinx x ii nx xx i f x ef x e C e e e=-=⋅==∑∑（）（，）. 14．用最小二乘法求解矛盾方程组2+314921x y x y x y =⎧⎪-=-⎨⎪-=-⎩. 解 方法一方程组可变形为31+22491122x y x y x y ⎧=⎪⎪-=-⎨⎪⎪-=-⎩，原问题转化成在已知三组离散数据3142211()922t f t ----下求一次最小二乘逼近函数1()P x x yt =+(x 与y 为一次函数的系数，t 为自变量)，取1H 基{}1,t ,求解法方程331133321113()()i i i i i i i i i i i x t f x t t t f x y =====⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦∑∑∑∑∑, 即3-3-93737-32x y ⎡⎤⎡⎤⎡⎤⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦，得到矛盾方程组的解为37=-3156=31x y ⎧⎪⎪⎨⎪⎪⎩. 方法二方程组可变形为31+22491122x y x y x y ⎧=⎪⎪-=-⎨⎪⎪-=-⎩，令(,)I x y 2223111=+-+4+9++2222x y x y x y --（）（）（）(,)I x y x ∂∂3111=2+-+24+9+2+2222x y x y x y ⨯⨯-⨯-（）（）（）=6618x y -+，(,)I x y y ∂∂331111=+44+9+222222x y x y x y ⨯--⨯--⨯-（）（）（） 37=3372x y -+- 令(,)0(,)0I x y x I x y y∂⎧=⎪∂⎪⎨∂⎪=⎪∂⎩， 得3373372x y x y -=-⎧⎪⎨-+-⎪⎩， 解之得矛盾方程组的解为37315631x y ⎧=-⎪⎪⎨⎪=⎪⎩. 习题47. 对列表函数 124810()0152127x f x求(5)(5).f f '''，解 一阶微商用两点公式(中点公式),得(8)(2)10(5),63f f f -'≈= 二阶微商用三点公式(中点公式),首先用插值法求(5)f , 由(4)5,(8)21,f f ==得一次插值函数1()411,L x x =-从而 1(5)(5)9f L ≈=,于是,2(2)2(5)(8)4(5).39f f f f -+''≈= 8. 导出数值数分公式)]23()2(3)2(3)23([1)(3)3(h x f h x f h x f h x f h x f ---++-+≈并给出余项级数展开的主部.解 由二阶微商的三点公式(中点公式),得213()[()2()()]2222h h h f x f x f x f x h h ''-≈+--+-,213()[()2()()]2222h h h hf x f x f x f x h ''+≈+-++-从而 (3)()()22()h h f x f x f x h''''+--≈3133=[()3()3()()]2222h h f x h f x f x f x h h +-++--- 将33()()()()2222h h f x h f x f x f x h ++--，，，分别在x 处展开,得2(3)3(4)4(5)55331313()=()()()()()()222!23!21313()()()()+()(1)4!25!2f x h f x f x h f x h f x h f x h f x h O h '''++⋅+⋅+⋅+⋅+⋅2(3)3(4)4(5)5511()=()()()()()()222!23!211()()()()()(2)4!25!2h h h h f x f x f x f x f x h h f x f x O h '''++⋅+⋅+⋅+⋅+⋅+2(3)3(4)4(5)5511()=()()()()()()()222!23!211()()()()()(3)4!25!2h h h h f x f x f x f x f x h h f x f x O h '''-+⋅-+⋅-+⋅-+⋅-+⋅-+2(3)3(4)4(5)55331313()=()()()()()()()222!23!21313()()()()()(4)4!25!2f x h f x f x h f x h f x h f x h f x h O h '''-+⋅-+⋅-+⋅-+⋅-+⋅-+(1)－(2)×3 +(3)×3－(4), 得(5)222131()[()2()()]()()22228h h h f x f x f x f x h f x h O h h ''--+--+-=-+,即余项主部为(5)21()8f x h -习 题 5 (P. 299)3. 设n n R A ⨯∈为对称矩阵，且011≠a ，经高斯消去法一步后，A约化为11120T a a A ⎡⎤⎢⎥⎣⎦，试证明2A 亦是对称矩阵. 证明设1111()=T ij aa A a A α⎛⎫= ⎪⎝⎭，其中 21311=n a a a α⎛⎫⎪ ⎪ ⎪ ⎪⎝⎭，121311=n a a a a ⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭，22232123=n n n nn a a a a a a a ⎛⎫⎪⎪ ⎪⎝⎭， 则经高斯消去法一步后，A 约化为111111110TT a a A a a α⎡⎤⎢⎥⎢⎥-⎢⎥⎣⎦， 因而211111T A A a a α=-，若n n R A ⨯∈为对称矩阵，则1A 为对称矩阵，且1=a α，易知211111T A A a a α=-为对称矩阵. 13. 设⎥⎦⎤⎢⎣⎡=989999100A(1) 计算2||||,||||A A ∞； (2) 计算∞)(A Cond ，及2)(A Cond . 解 (1) 计算||||=199A ∞,⎥⎦⎤⎢⎣⎡=989999100A，其特征值为1,299λ=，又⎥⎦⎤⎢⎣⎡=989999100A 为对称矩阵，则2=T A A A 的特征值为221,2(99λ=±，因此2||||99A ===+；(2)1989999100A --⎡⎤=-⎢⎥-⎣⎦，1||||=199A -∞， 所以1()=||||||||=9801Cond A A A -∞∞∞⋅，1989999100A --⎡⎤=-⎢⎥-⎣⎦为对称矩阵，其特征值为1,299λ=-± 则1112()=()T A A A ---的特征值为221,2(99λ=，因此12||||99A -===+所以1222()=||||||||Cond A A A -⋅2(99=+15. 设,n n n A R x R ⨯∈∈，求证 （1）1xx n x ∞∞≤≤； (2)∞∞≤≤An A A n11.证明 (2) 由（1）1x x n x∞∞≤≤，得1AxAx n Ax∞∞≤≤，则 11Ax Ax n Ax n x xx∞∞∞∞≤≤，从而11max max max nnnx Rx Rx RAxAx n Ax n xxx∞∞∀∈∀∈∀∈∞∞≤≤，由算子范数的定义max nx RAx Ax∞∞∀∈∞=，111max nx RAx A x∀∈=，得∞∞≤≤An A A n11.17. 设n n R W ⨯∈为非奇异阵，又设x为n R 上一向量范数，定义WxWx=，求证：Wx是nR 上向量的一种范数（称为向量的W 一范数）.证明 ①正定性，因Wx为一向量，0WxWx =≥，下证=0=0Wxx ⇔，⇒“”若=0Wx 即=0Wx ，由向量范数的正定性得=0Wx ，n n R W ⨯∈为非奇异阵，所以=0x ；⇐“”若=0x ，则=0Wx ，由向量范数的正定性得=0Wx 即=0Wx.②齐次性，任意实数α有=Wx W x Wxααα=，由向量范数的齐次性，得=WWxW x Wx Wx xααααα===；③ 三角不等式，任意实数,n n x R y R ∈∈，有+(+)=+Wx yW x y Wx Wy=，再由向量范数的三角不等式，得+(+)=+WWWx yW x y Wx Wy Wx Wy xy=≤+=+.习 题 6 (P.347)1. 设有方程组(b )1231231232211221x x x x x x x x x +-=⎧⎪++=⎨⎪++=⎩，考查用Jacobi迭代法,G-S 迭代法解此方程组的收敛性.解 系数矩阵分裂如下,122111221A -⎛⎫⎪= ⎪ ⎪⎝⎭D L U =--10022110112200-⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪=---- ⎪ ⎪⎪ ⎪ ⎪ ⎪--⎝⎭⎝⎭⎝⎭ Jacobi迭代矩阵为1()J D L U -=+=02211220-⎛⎫ ⎪-- ⎪ ⎪--⎝⎭， J 的特征方程为2211022λλλ-=，展开得 30λ=，即01λ=<，所以用Jacobi 迭代法解此方程组是收敛的.G-S 迭代矩阵为1()G D L U -=-11022=11012210--⎛⎫⎛⎫⎪ ⎪⋅- ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭100022=110010210-⎛⎫⎛⎫ ⎪ ⎪-⋅- ⎪⎪ ⎪ ⎪-⎝⎭⎝⎭122=023002-⎛⎫ ⎪- ⎪ ⎪⎝⎭， G 的特征方程为12221002λλλ---=-， 展开得 (1)(2)(2)0λλλ---=，即1λ=或2λ=，由迭代基本定理得用G-S 迭代法解此方程组是不收敛的.4. 设有方程组Ax b =，其中A 为对称正定阵，且有迭代公式(1)()()()k k k x x b Ax ω+=+- (0,1,k =)，试证明当20ωβ<<时，上述迭代法收敛(其中A 的特征值满足0()A αλβ<≤≤).证明 A 为对称正定阵， A 的特征值满足0()A αλβ<≤≤，且20ωβ<<，则0()2A ωλ<<又迭代公式可变形为(1)()()k k x I A x bωω+=-+ (0,1,k =)，从而迭代矩阵 B I A ω=-，迭代矩阵的特征值为1()A ωλ-，且满足11()1A ωλ-<-<，即 |()|1B λ<，由迭代基本定理得该迭代法是收敛的.5. 设111a a A aa a a⎛⎫ ⎪= ⎪ ⎪⎝⎭，其中a 为实数，试确定a 满足什么条件时，解Ax b =的Jacobi 迭代法收敛.解 系数矩阵分裂如下,111a a A aa a a ⎛⎫ ⎪= ⎪ ⎪⎝⎭D L U =--1001100a a aa aa--⎛⎫⎛⎫⎛⎫⎪ ⎪ ⎪=---- ⎪ ⎪⎪ ⎪ ⎪ ⎪--⎝⎭⎝⎭⎝⎭Jacobi迭代矩阵为1()J D L U -=+=000a a aa a a--⎛⎫⎪-- ⎪ ⎪--⎝⎭，J 的特征方程为0a aa a aaλλλ=，展开得 323320a a λλ--=，即a λ=-或2a λ=-，()max{||,|2|}J a a ρ=--()1J ρ<当且仅当1122a -<<，所以当1122a -<<时，解Ax b=的Jacobi 迭代法收敛.。
数值分析上机实习题及答案.docx
方詡文金兴:爭［数值分析］2017-2018第二学期上机实习题1：编程计算亍丄，其中C= 4. 4942x10307,给出并观察计算结心C"果，若有问题，分析之。
解:mat lab 编程如下:E) funct ion diy i ti formatlong g;n 二input C 输入ii 值= c= 4.4942E307; sum 二0; s 二 0;E3 for i = l:n s = l/ (c*i);>> diyiti 输入n 值:10 104.6356e-308 >> diyiti输入ri 值:1001004.6356e-308 >> diyiti 输入n 值:1000 10004.6356e-308 >> diyiti揄入n 值* 1000001000004・ 6356e-308 >> diyiti输入n 值；1000000001000000004.6356e-308图二：运行结果Mat lab 中，forma t long g 对双精度，显示15位定点或浮点格式，由上图 可知，当输入较小的n 值5分别取10, 100, 1000, 100000, 100000000）的时候, 结果后面的指数中总是含有e-308,这和题目中的C 值很相似，我认为是由于分 母中的C 值相对于n 值过大，出现了 “大数吃小数”的现彖，这是不符合算法原 则的。
2：利用牛顿法求方程-1^ = 2于区间241的根，考虑不同初值下牛顿法的收敛情况。
解：牛顿法公式为:利用mat lab 编程function di2ti21 3i=l ；2 2.85208156699784 xO 二input （'输入初值x0：‘ ）;A 二[i x0];3 2.55030468822809 t=x0+ (x0-log (xO) -2) /(1-1/xO) ; %迭代函数4 1.91547247100476 三 while (abs (t _x0)>0.01)i=i+l; 5 0.37867158538991 xO 二 t; 6 0.774964959780275 A = [A;i xO];t =x0+(x0-log(xO)-2)/(1-1/xO): 7 4.11574081641933 cnd| 8 5.04162436446126 disp (A);96.81782826645596当输入初值二3的时候并不能收敛。
数值分析课后部分习题答案
解
x * = 2.00021 = 0.200021 × 101 ，即 m = 1
1 1 × 10m − n = × 10−3 ， 2 2
由有效数字与绝对误差的关系得 即
m − n = −3 ，所以， n = 2 ； y* = 0.032 = 0.32 × 101 ，即 m = 1
由有效数字与绝对误差的关系得 即
m − n = −3 ，所以， n = 4 ； z * = 0.00052 = 0.52 × 10−3 ，即 m = −3
1 1 × 10m − n = × 10−3 ， 2 2
由有效数字与绝对误差的关系得 即
m − n = −3 ，所以， n = 0 .
1 1 × 10m − n = × 10−3 ，Fra bibliotek2 2=
f [x1 , x2 ,⋯ , x n ]-f [ x0 , x1 ,⋯ , x n−1 ] g[ x1 , x2 ,⋯ , x n ] − g[ x0 , x1 ,⋯ , x n−1 ] + x n − x0 x n − x0
( x − 1)( x − 2)( x − 3) 1 =- ( x − 1)( x − 2)( x − 3) ， (0 − 1)(0 − 2)(0 − 3) 6
x ( x − 2)( x − 3) 1 = x ( x − 2)( x − 3) ， (1 − 0)(1 − 2)(1 − 3) 2 x( x − 1)( x − 3) 1 =- x( x − 1)( x − 3) ， (2 − 0)(2 − 1)(2 − 3) 2 x( x − 1)( x − 2) 1 = x ( x − 1)( x − 2) ， (3 − 0)(3 − 1)(3 − 2) 6
研究生数值分析答案_2
● 当22()3x x ϕ+=时，'12()3x x ϕ=，因此'1(2) 1.3333ϕ=>1，，因此，该迭代格式不收敛。
● 当2()x ϕ='2()x ϕ=，因此'2(2)0.75ϕ=<1，，因此，该迭代格式收敛。
● 当32()3x x ϕ=-时，'322()x xϕ=，因此'3(2)0.5ϕ=<1，，因此，该迭代格式收敛。
● 当242()23x x x ϕ-=- 时，2'44224()1213x x x x x ϕ--=-+，因此'4(2)0ϕ=<1，，因此，该迭代格式收敛。
(2)、● 当22()3x x ϕ+=时，迭代法计算公式是20122.5,3k k x x x ++==，程序如下： >> fi=inline('(x.*x+2)/3');x0=2.5;er=1;k=0;while er>0.00001x=fi(x0);er=abs(x-x0);x0=xk=k+1end运行结果如下：x0 =2.7500k =1x0 =3.1875k =2x0 =4.0534k =3x0 =6.1433k =x0 =13.2468k =5x0 =59.1589k =6x0 =1.1673e+003 k =7x0 =4.5416e+005 k =8x0 =6.8755e+010 k =9x0 =1.5757e+021 k =10x0 =8.2765e+041 k =11x0 =2.2834e+083 k =12x0 =1.7379e+166 k =13x0 =Infk =14x0 =Infk =15由以上计算结果看，序列是发散的，运行14次已经超出计算机的识别范围，当2()x ϕ迭代法计算公式是1k x +=程序运行结果如下：>> fi=inline('sqrt(3*x-2)');x0=2.5;er=1;k=0;while er>0.00001x=fi(x0);er=abs(x-x0);x0=xk=k+1;endx0 =2.3452x0 =2.2440x0 =2.1753x0 =2.1274x0 =2.0934x0 =2.0689x0 =2.0510x0 =2.0379x0 =2.0282x0 =2.0211x0 =2.0157x0 =2.0118x0 =2.0088x0 =2.0066x0 =2.0049x0 =2.0037x0 =2.0028x0 =2.0021x0 =2.0016x0 =2.0012x0 =2.0009x0 =2.0007x0 =2.0005x0 =2.0004x0 =2.0003x0 =2.0002x0 =2.0002x0 =2.0001x0 =2.0001x0 =2.0001x0 =2.0000x0 =2.0000x0 =2.0000>>由以上计算结果看，序列收敛与2，所以x=2是f(x)= 232x x -+=0的根。
数值分析课后习题答案
x2 6.6667x2 8.205
再解
1
15 56
x31.785,7得 x35.769
1 25069x4 0.47847x4 1.4872
1 x5 5.3718 x5 5.3718
2-10.证明下列不等式:
(1)x-yx-z+z-y; (2)|x-y|x-y;
证明 (1)x-y=(x-z)+(z-y)x-z+z-y
b.用Gauss消元法
102 x y 1 x y 2
回代得解: y=1, x=0.
102 x Байду номын сангаасy 1
100y 100
再用列主元Gauss消元法
102 x y 1 x y 2
回代得解: y=1, x=1.
x y
y 1
2
2-8.用追赶法求解方程组:
4 1
x1 100
1 4 1
x2 0
3-8.判定求解下列方程组的SOR方法的收敛性.
2 1 0 0 x1 1
1
0 0
2 1 0
1 2 1
0 12
x2 x3 x4
0 00
解 直接可验证系数矩阵A是负定矩阵,所以-A是对称
1-3.为了使101/2的相对误差小于0.01%,试问应取几位 有效数字?
解 因为101/2=3.162…=0.3162…10,若具有n位有效 数字,则其绝对误差限为0.5 101-n ,于是有
r=0.5101-n/3.162…<0.5101-n/3<0.01% 因此只需n=5.即取101/2=3.1623
1 2
0
12 1,
1 2
1 2
0
12
数值分析上机题参考答案.docx
如有帮助欢迎下载支持数值分析上机题姓名：陈作添学号： 040816习题 120．（上机题）舍入误差与有效数N11 3 1 1设S N，其精确值为 。
22 2 NN 1j 2j1（1）编制按从大到小的顺序111 ，计算 S 的通用程序。
S N1 321N 21 N22（2）编制按从小到大的顺序111，计算 S 的通用程序。
S N1(N 1)2 122 1NN 2（3）按两种顺序分别计算S 102 ， S 104 ， S 106 ，并指出有效位数。
（编制程序时用单精度）（4）通过本上机题，你明白了什么？按从大到小的顺序计算 S N 的通用程序为： 按从小到大的顺序计算 S N 的通用程序为：#include<iostream.h> #include<iostream.h> float sum(float N) float sum(float N) {{float j,s,sum=0; float j,s,sum=0; for(j=2;j<=N;j++) for(j=N;j>=2;j--) {{s=1/(j*j-1); s=1/(j*j-1); sum+=s;sum+=s;}}return sum;return sum;}}从大到小的顺序的值从小到大的顺序的值精确值有效位数从大到小从小到大0.7400490.740050.74004965 S 1020.7498520.74990.74994 4S 1040.7498520.7499990.74999936S 106通过本上机题， 看出按两种不同的顺序计算的结果是不相同的，按从大到小的顺序计算的值与精确值有较大的误差， 而按从小到大的顺序计算的值与精确值吻合。
从大到小的顺序计算得到的结果的有效位数少。
计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低，我们在计算机中进行同号数的加法时，采用绝对值较小者先加的算法，其结果的相对误差较小。
数值分析课后部分习题答案
证明 由差商的定义 (a) 如果 F ( x ) = cf ( x ) ，则
F [ x0 , x1 ,⋯ , xn ] =
=
F [ x1 , x2 ,⋯ , xn ]-F [ x0 , x1 ,⋯ , xn− 1 ] x n − x0
cf [ x1 , x2 , ⋯ , xn ]-cf [ x0 , x1 ,⋯ , xn −1 ] x n − x0 f [ x1 , x2 , ⋯ , xn ]-f [ x0 , x1 ,⋯ , xn−1 ] = cf [ x0 , x1 , ⋯ , xn ] . x n − x0
1 1 1 1 |e( x*)| ≤ × 10m − n = × 10−2 ， |e( y*)| ≤ × 10m − n = × 10 −2 ， 2 2 2 2 1 1 |e( z*)| ≤ × 10 m − n = × 10 −2 ， 2 2 | e( y * z*) |≈| z * e ( y*) + y * e ( z *) |≤ z * | e ( y *) | + y * | e (z *) |
m − n = −3 ，所以， n = 4 ； z * = 0.00052 = 0.52 × 10−3 ，即 m = −3
1 1 × 10m − n = × 10−3 ， 2 2
由有效数字与绝对误差的关系得 即
m − n = −3 ，所以， n = 0 .
1 1 × 10m − n = × 10−3 ， 2 2
1 1 ≤ 2.35 × × 10−2 + 1.84 × × 10−2 = 2.095 × 10−2 ， 2 2 1 | e( x * + y * z*) |≈| e( x*) + e( y * z*) |≤ × 10 −2 + 2.095 × 10−2 2 1 = 0.2595 × 10−1 ≤ × 10−1 ， 2
数值分析答案Ch1-Ch2
第一二章习题提示第一章误差1.计算球体积要使相对误差限为1%，问度量半径R 时允许的相对误差是多少？答案：V =43πR 3，εr (V )=ε(V )V =4πR 2ε(R )43πR 3=3εr (R )。
球体积要使相对误差限为1%，R 时允许的相对误差是1300=0.003333。
2.考虑正弦函数sin x 的求值，特别是数据传递误差，即自变量x 发生扰动h 时函数值的误差。
(1)估计sin x 的绝对误差。
(2)估计sin x 的相对误差。
(3)估计这个问题的条件数。
(4)自变量x 为何值时，这个问题高度敏感？答案：(1)cos xh 。
(2)cot xh 。
(3)相对条件数x cot x 。
(4)x =kπ，k =0时，这个问题高度敏感。
注：误差估计与近似值取法有关，要舍去高阶无穷小。
如：sin (x +h )≈sin x +cos xh ，则误差为−12sin xh 2等。
3.设Y 0=28，按递推公式Y n =Y n −1−1100√783，(n =1,2,...)计算到Y 100。
若取√783≈27.982(保留5位有效数字)，试问计算Y 100将有多大误差。
答案：Y 100=Y 0−√783，ε(Y 100)=ε(27.982)(假设Y 0无误差)。
Y 100的误差限为12×10−3。
4.正方形的边长大约为100cm ，问测量时允许多大的误差才能使其面积误差不超过1cm 2。
答案：测量边长的误差应不超过0005cm 。
上机题1.编程观察无穷级数∞∑n =11n的求和计算。
(1)采用IEEE 单精度浮点数，观察当n 为何值时求和结果不再变化，将它与理论分析的结论进行比较(注：在MATLAB 中可用single 命令将变量转成单精度浮点数)。
(2)用IEEE 双精度浮点数计算(1)中前n 项的和，评估IEEE 单精度浮点数计算结果的误差。
答案：format longs =single(0);u =single(1);t =0;k =1;while u ~=su =s;s =s +1/k;t =t +1/k;k =k+1;endk -1s -t这里单精度浮点数的和为s =15.404，双精度浮点数的和为t =15.133。
数值分析上机试题对应参考答案
一、 问答题1、什么是近似值x * 有效数字？若近似值x*的误差限是某一位的半个单位，该位到x*的第一位非零数字共有n 位，就说x*有n 位有效数字。
它可表示为X=±10m ×(a 1+a 2×10-1+…+a n ×10-(n-1),其中a i (i=1,2,…,n)是0到9中的一个数字，a 1≠0,m 为整数，且︱x －x *︱≠21×10m-n+12、数值计算应该避免采用不稳定的算法，防止有效数字的损失. 因此，在进行 数值运算算法设计过程中主要注意什么？ （1）简化计算过程，减少运算次数； （2）避免两个相近的数相减；（3）避免除数的绝对值远小于被除数的绝对值； （4）防止大数“吃掉”小数的现象；（5）使用数值稳定的算法，设法控制误差的传播。
3、写出“n 阶阵A 具有n 个不相等的特征值”的等价条件（至少写3 个）(1)|A|不为零(2)n 阶矩阵A 的列或行向量组线性无关 (3)矩阵A 为满秩矩阵(4)n 阶矩阵A 与n 阶可逆矩阵B 等价4、迭代法的基本思想是什么？就是用某种极限过程去逐步逼近线性方程组精确解得方法。
其基本思想为：先任取一组近似解初值X 0，然后按照某种迭代原则，由X 0计算新的近似解X 1，以此类推，可计算出X 2,X 3，…X K ，。
,如果｛X ｝收敛，则取为原方程组的解。
5、病态线性方程组的主要判断方法有哪些？（1）系数矩阵的某两行（列）几乎近似相关 （2）系数矩阵的行列式的值很小（3）用主元消去法解线性方程组时出现小主元（4）近似解x*已使残差向量r=b-Ax*的范数很小，但该近似解仍不符合问题要求。
6、Lagrange 插值的前提条件是什么？并写出二次Lagrange 插值的基函数。
前提条件是：⎩⎨⎧≠==i j i j x j，，（01)l i .2,1,0,n j i ， = 二次Lagrange 插值的基函数：()))(())((2010210x x x x x x x x x l ----=()))(())((2101201x x x x x x x x x l ----= ()))(())((1202102x x x x x x x x x l ----=7、什么是数值积分的代数精度？如果某一个求积公式对于次数不超过m 的多项式均能准确地成立，但对于m+1次多项式就不准确成立，则称该求积公式具有m 次代数精度（或代数精确度）。
数值分析课后题答案
数值分析 第二章2．当1,1,2x =-时，()0,3,4f x =-,求()f x 的二次插值多项式。
解：0120121200102021101201220211,1,2,()0,()3,()4;()()1()(1)(2)()()2()()1()(1)(2)()()6()()1()(1)(1)()()3x x x f x f x f x x x x x l x x x x x x x x x x x l x x x x x x x x x x x l x x x x x x x ==-===-=--==-+-----==------==-+--则二次拉格朗日插值多项式为220()()k k k L x y l x ==∑0223()4()14(1)(2)(1)(1)23537623l x l x x x x x x x =-+=---+-+=+- 6．设,0,1,,j x j n =为互异节点，求证：（1）0()nkkj j j x l x x=≡∑ (0,1,,);k n =（2）0()()0nk jj j xx l x =-≡∑ (0,1,,);k n =证明（1） 令()kf x x = 若插值节点为,0,1,,j x j n =，则函数()f x 的n 次插值多项式为0()()nk n j j j L x x l x ==∑。
插值余项为(1)1()()()()()(1)!n n n n f R x f x L x x n ξω++=-=+又,k n ≤(1)()0()0n n f R x ξ+∴=∴=0()nk kj j j x l x x =∴=∑ (0,1,,);k n =000(2)()()(())()()(())nk j j j n nj i k i k j j j i nnik ii kj j i j x x l x C x x l x C x x l x =-==-==-=-=-∑∑∑∑∑0i n ≤≤又 由上题结论可知()nk ij jj x l x x ==∑()()0ni k i ik i k C x x x x -=∴=-=-=∑原式∴得证。
数值分析课后习题集部分参考答案解析
数值分析课后习题部分参考答案Chapter 1（P10）5. 求2的近似值*x ，使其相对误差不超过%1.0。
解： 4.12=。
设*x 有n 位有效数字，则nx e -⨯⨯≤10105.0|)(|*。
从而，1105.0|)(|1*nr x e -⨯≤。
故，若%1.0105.01≤⨯-n，则满足要求。
解之得，4≥n 。
414.1*=x 。
（P10）7. 正方形的边长约cm 100，问测量边长时误差应多大，才能保证面积的误差不超过12cm 。
解：设边长为a ，则cm a 100≈。
设测量边长时的绝对误差为e ，由误差在数值计算的传播，这时得到的面积的绝对误差有如下估计:e ⨯⨯≈1002。
按测量要求，1|1002|≤⨯⨯e 解得，2105.0||-⨯≤e 。
Chapter 2（P47）5. 用三角分解法求下列矩阵的逆矩阵：⎪⎪⎪⎭⎫ ⎝⎛--=011012111A 。
解：设()γβα=-1A。
分别求如下线性方程组：⎪⎪⎪⎭⎫ ⎝⎛=001αA ，⎪⎪⎪⎭⎫ ⎝⎛=010βA ，⎪⎪⎪⎭⎫⎝⎛=100γA 。
先求A 的LU 分解（利用分解的紧凑格式），⎪⎪⎪⎭⎫⎝⎛-----3)0(2)1(1)1(2)0(1)1(2)2(1)1(1)1(1)1(。
即，⎪⎪⎪⎭⎫ ⎝⎛=121012001L ，⎪⎪⎪⎭⎫⎝⎛---=300210111U 。
经直接三角分解法的回代程，分别求解方程组，⎪⎪⎪⎭⎫ ⎝⎛=001Ly 和y U =α，得，⎪⎪⎪⎭⎫ ⎝⎛-=100α；⎪⎪⎪⎭⎫ ⎝⎛=010Ly 和y U =β，得，⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=323131β；⎪⎪⎪⎭⎫ ⎝⎛=100Ly 和y U =γ，得，；⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=313231γ。
所以，⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛---=-3132132310313101A 。
（P47）6. 分别用平方根法和改进平方根法求解方程组：⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----816211515311401505231214321x x x x 解: 平方根法：先求系数矩阵A 的Cholesky 分解（利用分解的紧凑格式），⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----1)15(2)1(1)5(3)3(3)14(2)0(1)1(1)5(2)2(1)1(，即，⎪⎪⎪⎪⎪⎭⎫⎝⎛--=121332100120001L ，其中，TL L A ⨯=。
数值分析上机题目详解
第一章一、题目设∑=-=NN j S 2j 211，其精确值为)11123(21+--N N 。
1) 编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算S N 的通用程序。
2) 编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算S N 的通用程序。
3) 按两种顺序分别计算64210,10,10S S S ,并指出有效位数。
（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？二、通用程序N=input('Please Input an N (N>1):'); AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); Sn1=single(0);for a=2:N; Sn1=Sn1+1/(a^2-1); endSn2=single(0);for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); endfprintf('The value of Sn (N=%d)\n',N);fprintf('Accurate Calculation %f\n',AccurateValue); fprintf('Caculate from large to small %f\n',Sn1); fprintf('Caculate from small to large %f\n',Sn2); disp('____________________________________________________')三、结果从结果可以看出有效位数是6位。
感想：可以得出，算法对误差的传播有一定的影响，在计算时选一种好的算法可以使结果更为精确。
从以上的结果可以看到从大到小的顺序导致大数吃小数的现象，容易产生较大的误差，求和运算从小数到大数所得到的结果才比较准确。
数值分析 第三章上机答案
1、程序：x=-1:0.2:1;f=1./(1+25*x.^2);y=polyval(f,x); %计算原函数每个x 值所对应的函数值 p=polyfit(x,y,3) %对（x ，y ）进行三次拟合并输出三次多项式系数 z=polyval(p,x); %计算拟合后的多项式对应x 点的函数值 plot(x,y,'r',x,z,'b') %画图输出：p =1.8222 1.2090 -0.3619 -0.1140则三次曲线拟合的方程为：1140.03619.02090.18222.1)(233--+=x x x x p2、3次和4次多项式拟合：程序：x=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];y=[1.0 0.41 0.50 0.61 0.91 2.02 2.46];%3次拟合%p3=polyfit(x,y,3)xi=0:0.1:1.0;yi=polyval(p3,xi);subplot(1,2,1);plot(x,y,'*',xi,yi,'r');xlabel('x');ylabel('y');%4次拟合%p4=polyfit(x,y,4)xi=0:0.1:1.0;yi=polyval(p4,xi);subplot(1,2,2);plot(x,y,'*',xi,yi,'r');xlabel('x');ylabel('y');输出：p3 =-6.6221 12.8147 -4.6591 0.9266p4 =2.8853 -12.3348 16.2747 -5.2987 0.9427则3次拟合多项式为：9266.06591.48147.126221.6)(233+-+-=x x x x p 4次拟合多项式为：9427.02987.52747.163348.128853.2)(2344+-+-=x x x x x p另一函数拟合：定义函数：function [C,L]=lagran(x,y)%x 为n 个节点的横坐标所组成的向量，y 为纵坐标所组成的向 %C 为所求的牛顿插值多项式的系数构成w=length(x);n=w-1;L=zeros(w,w);for k=1:n+1V=1;for j=1:n+1if k~=jV=conv(V,poly(x(j)))/(x(k)-x(j));% conv 求积，poly(x)将该多项式的系数赋给向量endendL(k,:)=V;endC=y*L输入命令：x=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];y=[1.0 0.41 0.50 0.61 0.91 2.02 2.46];cc=polyfit(x,y,4);xx=x(1):0.1:x(length(x));yy=polyval(cc,xx);plot(xx,yy,'r');hold on ;plot(x,y,'x');xlabel('x');ylabel('y');x=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];y=[0.94 0.58 0.47 0.52 1.00 2.00 2.46];%y 中的值是根据上面两种拟合曲线而得到的估计数据，不是真实数据[C,L]=lagran(x,y);xx=0:0.01:1.0;yy=polyval(C,xx);hold onplot(xx,yy,'b',x,y,'.');输出：C =40.6746 -110.2183 110.3671 -57.3264 23.4994 -5.4764 0.9400。
数值分析实验(参考答案)
数值分析上机实验学生姓名： 学号： 教师：实验1：1. 实验项目的性质和任务通过上机实验，使学生对病态问题、线性方程组求解和函数的数值逼近方法有一个初步理解。
2．教学内容和要求 1）对高阶多多项式201()(1)(2)(20)()k p x x x x x k ==---=-∏编程求下面方程的解 19()0p x x ε+=并绘图演示方程的解与扰动量ε的关系。
Matlab 程序：x=1:20;y=zeros(1,20); ve=zeros(1,21); plot(x,y,'o') hold on ; pause; for x=1:5ve(2)=10^(-x);e=roots(poly(1:20)+ve);plot(e,'*') hold on pause; end0510152025303540-20-15-10-5051015202）对2~20n =，生成对应的Hilbert 矩阵，计算矩阵的条件数；通过先确定解获得常向量b 的方法,确定方程组 n H x b =最后，用矩阵分解方法求解方程组，并分析计算结果。
Matlab 程序：for n=2:20; H=hilb(n); ca=cond(H,2) x=(1:n)'; b=H*x; [L,U]=lu(H); y=L\b;x1=U\yplot(x,x,'o',x1,x1,'*') hold on pause; end-500-400-300-200-100100200300400500-500-400-300-200-10001002003004005003）对函数 21()[1,1]125f x x x=∈-+的Chebyshev 点 (21)cos()1,2,...,12(1)kk x k n n π-==++编程进行Lagrange 插值，并分析插值结果。
Matlab 程序：function y=lagrangen(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i);s=0;for k=1:nL=1;for j=1:nif j~=kL=L*(z-x0(j))/(x0(k)-x0(j));endends=s+L*y0(k);endy(i)=s;endy;for n=5:20x=-1:0.01:1; y=1./(1+25*x.^2);plot(x,y,'r')hold onk=n+1:-1:1;x0=cos((2*k-1)*pi./(2*(n+1))),y0=1./(1+25*x0.^2);x=-1:0.01:1; y1=lagrangen(x0,y0,x);plot(x0,y0,'o',x,y1), pausehold off end-1-0.8-0.6-0.4-0.200.20.40.60.81-0.200.20.40.60.811.2-1-0.8-0.6-0.4-0.200.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81。
数值分析课后习题部分参考答案.
数值分析课后习题部分参考答案Chapter 1（P10）5. 求2的近似值*x ，使其相对误差不超过%1.0。
解： 4.12=。
设*x 有n 位有效数字，则n x e -⨯⨯≤10105.0|)(|*。
从而，1105.0|)(|1*nr x e -⨯≤。
故，若%1.0105.01≤⨯-n，则满足要求。
解之得，4≥n 。
414.1*=x 。
（P10）7. 正方形的边长约cm 100，问测量边长时误差应多大，才能保证面积的误差不超过12cm 。
解：设边长为a ，则cm a 100≈。
设测量边长时的绝对误差为e ，由误差在数值计算的传播，这时得到的面积的绝对误差有如下估计:e ⨯⨯≈1002。
按测量要求，1|1002|≤⨯⨯e 解得，2105.0||-⨯≤e 。
Chapter 2（P47）5. 用三角分解法求下列矩阵的逆矩阵：⎪⎪⎪⎭⎫ ⎝⎛--=011012111A 。
解：设()γβα=-1A。
分别求如下线性方程组：⎪⎪⎪⎭⎫ ⎝⎛=001αA ，⎪⎪⎪⎭⎫ ⎝⎛=010βA ，⎪⎪⎪⎭⎫ ⎝⎛=100γA 。
先求A 的LU 分解（利用分解的紧凑格式），⎪⎪⎪⎭⎫ ⎝⎛-----3)0(2)1(1)1(2)0(1)1(2)2(1)1(1)1(1)1(。
即，⎪⎪⎪⎭⎫ ⎝⎛=121012001L ，⎪⎪⎪⎭⎫⎝⎛---=300210111U 。
经直接三角分解法的回代程，分别求解方程组，⎪⎪⎪⎭⎫ ⎝⎛=001Ly 和y U =α，得，⎪⎪⎪⎭⎫ ⎝⎛-=100α；⎪⎪⎪⎭⎫ ⎝⎛=010Ly 和y U =β，得，⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=323131β；⎪⎪⎪⎭⎫ ⎝⎛=100Ly 和y U =γ，得，；⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=313231γ。
所以，⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛---=-3132132310313101A。
（P47）6. 分别用平方根法和改进平方根法求解方程组：⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----816211515311401505231214321x x x x 解:平方根法：先求系数矩阵A 的Cholesky 分解（利用分解的紧凑格式），⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----1)15(2)1(1)5(3)3(3)14(2)0(1)1(1)5(2)2(1)1(，即，⎪⎪⎪⎪⎪⎭⎫⎝⎛--=121332100120001L ，其中，TL L A ⨯=。
数值分析上机答案
第六章 2、证明：对于矩阵A 进行分解有：000300002000020001210002000A L D U=++-⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪=++ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭那么：对于Jacobi 迭代其迭代矩阵为：11()20033000021020001020022101102J B D L U --=+⎛⎫- ⎪-⎛⎫⎛⎫ ⎪ ⎪ ⎪⎪=+= ⎪ ⎪⎪ ⎪ ⎪- ⎪⎝⎭⎝⎭⎪- ⎪⎝⎭对于Gauss-Seidel 迭代其迭代矩阵为：11()20033000021020001002212000110012g s J B B D L U---==-⎛⎫- ⎪-⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪=+= ⎪ ⎪ ⎪ ⎪ ⎪- ⎪⎝⎭⎝⎭ ⎪ ⎪⎝⎭（1）对于B J 我们求()J B ρ有：223111det(I )0()212112J B λλλλλλ-=-=--所以12,30,λλ==则有：()1J B ρ=< ，所以Jacobi 迭代收敛。
对于Gauss-Seidel 迭代，求其()g s B ρ-有：11()112g s g sB B ρ--∞<=< 所以Gauss-seidel 迭代收敛。
（2）并且有11()()12g s J B B ρρ-=<= 所以Gauss-Seidel 迭代收敛速度更快。
3、解：求取系数矩阵并进行分解有：4101143,50143A b -⎛⎫⎛⎫ ⎪ ⎪=--= ⎪ ⎪ ⎪ ⎪--⎝⎭⎝⎭那么SOR 迭代公式为：取初始值为（0,0,0）有：(0)1(1)()(1)()1(0,0,0)()/1,2,3T i n k k k k ii i ij i ij i ii j j i x x x b a x a x a i ω-++==⎧=⎪⎪=+--⎨⎪⎪=⎩∑∑ 那么有：对于ω=1.03有：3(1)(0)(0)111113(1)(0)(1)(0)22221122222(1)(0)(1)(0)33333333311.03()/40.25751.03()/ 1.353806251.03()/0.423894890625j j j j j j j j j x x b a x x x b a x a x a x x b a x a x a ===⎧=+-=⎪⎪⎪⎪=+--=⎨⎪⎪=+--=-⎪⎪⎩∑∑∑ 那么*(1)0.35380625x x ∞-=，条件不满足，继续迭代有：6次。
数值分析第五版上机实验答案实验一~实验六
实验一Lagrange插值算法实验目的：掌握拉格朗日（Lagrange）插值算法的基本原理，理解插值基函数的性质，掌握基本的误差概念。
学习用计算机语言编写程序实现算法。
[参考程序]#include "stdio.h"//定义插值节点及所求点数据，根据题目不同而修改double x[]={0.32,0.34,0.36};double y[]={0.314567,0.333487,0.352274};double xx=0.3367;// Lagrange插值算法函数，利用循环计算具有对称性的基函数和最终结果double Lagrange(double xxx,int n){int i;double result=0,temp;for(i=0;i<=n;i++){temp=1;for(int j=0;j<= n;j++){if(j!=i){temp=temp*(xxx-x[j])/(x[i]-x[j]);}}result=result+temp*y[i];}return result;}void main(){int n;printf("Please input n:");scanf("%d",&n);printf("Sin(%f) = %f \n",xx,Lagrange(xx,n));}实验二Newton均差插值算法实验目的：掌握Newton均差插值算法的基本原理，理解均差的概念，掌握均差表的计算方法。
学习用计算机语言编写程序实现算法。
[参考程序]#include "stdio.h"#define N 10double f[N][N];//定义插值节点及所求点数据，根据题目不同而修改double x[]={0.4,0.55,0.65,0.80,0.90,1.05};double y[]={0.41075,0.57815,0.69675,0.88811,1.02652,1.25382};double fx(int i,int j);double S(int start,int end,double xx);main(){int loopi,loopj,n;double result,xx;scanf("%d",&n);scanf("%lf",&xx);for(loopi=0;loopi<=n;loopi++){//零阶均差作为均差表二维数组的第0列f[loopi][0]=y[loopi];}//循环计算均差表中的所有数据for(loopi=1;loopi<=n;loopi++){for(loopj=1;loopj<=loopi;loopj++){f[loopi][loopj]=fx(loopi,loopj);}}result=S(0,n,xx);printf("Result is: %f",result);return 1;}//求均差的函数double fx(int i,int j){if(j==0){return f[i][j];}else{//这种表示方法需要注意两个x的下标return (fx(i,j-1)-fx(i-1,j-1))/(x[i]-x[i-j]);}}//用秦九韶算法计算插值多项式结果double S(int start,int end,double xx){if(start==end){return f[end][end];}else{return (S(start+1,end,xx)*(xx-x[start])+f[start][start]);}}实验三Newton差分插值算法实验目的：掌握Newton差分插值算法的基本原理，理解差分的概念，掌握差分表的计算方法。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2009级研究生《数值分析》上机作业院系电气工程学院专业控制理论与控制工程姓名马凯指导教师代新敏2009年12月29日第一题（二问）：超松弛法求方程组根1.解题理论依据或方法应用条件:超松弛算法是在GS 方法已求出x (m),x (m-1)的基础上，经过重新组合得到新序列。
如能恰当选择松弛因子ω，收敛速度会比较快。
当ω>1时，称为超松弛法，可以用来加速收敛。
其具体算法为：)()1(1)1(11)()1()(i ni j m jij i j m j ij m i m i g xb x b x x +++-=∑∑+=--=-ωω2.计算程序（使用软件：VC ）： #include<stdio.h> #define w 1.4 main(){float a[10][10]={ {0,0,0,0,0,0,0,0,0,0},{0,12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743}, {0,2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124}, {0,-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103}, {0,1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585}, {0,-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317}, 0,0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417}, {0,1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474},{0,3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782}, {0,-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001}}; float b[10][1]={{0},{2.1874369},{33.992318},{-25.173417},{0.84671695},{1.784317},{-86.612343},{1.1101230},{4.719345},{-5.6784392}};float x[10][10]={{0},{0},{0},{0},{0},{0},{0},{0},{0},{0}}; /*由x(0)=0得到其第一列全为零*/ float sum1=0,sum2=0; int i,m,j;for(m=1;m<=9;m++) for(i=1;i<=9;i++) {sum1=0; for(j=1;j<=(i-1);j++)sum1+=(-a[i][j]/a[i][i])*x[j][m]; /*计算第一个累加和*/ sum2=0; for(j=(i+1);j<=9;j++)sum2+=(-a[i][j]/a[i][i])*x[j][m-1]; /*计算第二个累加和*/ x[i][m]=(1-w)*x[i][m-1]+w*(sum1+sum2+b[i][0]/a[i][i]); /*用SOR 方法计算*/ }printf("x1为：%lf\n",x[1][9]); printf("x2为：%lf\n",x[2][9]); printf("x3为：%lf\n",x[3][9]); printf("x4为：%lf\n",x[4][9]); printf("x5为：%lf\n",x[5][9]); printf("x6为：%lf\n",x[6][9]);printf("x7为：%lf\n",x[7][9]); printf("x8为：%lf\n",x[8][9]); printf("x9为：%lf\n",x[9][9]); }3.计算结果4.问题讨论（误差分析、上机出现情况等）这道题目是所有题目中编写最顺利的，一次即顺利得出结果，当然这道题目还是有应该注意到地方，一是注意两个求和的清零，二是注意下标，不要弄混行标和列标。
第一题（三问）：列主元素法求方程组根1.解题理论依据或方法应用条件:所谓列主元消去法是，对矩阵作恰当的调整，选取绝对值最大的元素作为主元素。
然后把矩阵化为上三角阵，再进行回代，求出方程的解。
算法为： ][]1[1i b i d =+ , ]][1[]2[1i i B i a +=+ ,]][[]1[1i i B i b =+ , ]1][[]1[1+=+i i B i c , q1[0]=0 , u1[0]=0 ,8,,2,1]),[1][1][1(][1][1⋅⋅⋅=⋅+-=i i q i a i b i c i q9,,2,1]),1[1][1][1(])[1][1][1(}[1⋅⋅⋅=-⋅+⋅-=i i q i a i b i u i a i d i u x[9]=u1[9]1,,7,8],[1]1[][1][⋅⋅⋅=++⋅=i i u i x i q i x2.计算程序（使用软件：VC）：#include<stdio.h>#include<math.h>main(){double a[10][11]={ {0,0,0,0,0,0,0,0,0,0,0},{0,12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743,2. 1874369},{0,2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124,3 3.992318},{0,-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,-25.173417},{0,1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585,0. 84671695},{0,-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317,1. 784317},{0,0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417,-8 6.612343},{0,1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.21347 4,1.1101230},{0,3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782,4.7 19345},{0,-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001,-5 .6784392}};double x[10],y[11];double sum=0,max;int i,j,j1,m,n,k,g,h;for(j=1;j<=8;j++){max=a[j][j];for(k=j+1;k<=9;k++){if(fabs(a[k][j])>fabs(max))max=a[k][j];}for(k=j;k<=9;k++){if(a[k][j]/max==1)g=k;} /*确定本列绝对值最大元素所在行：g行*/ for(h=j;h<=10;h++)y[h]=a[j][h];for(h=j;h<=10;h++)a[j][h]=a[g][h];for(h=j;h<=10;h++)a[g][h]=y[h]; /*将g行与J行所有元素进行交换*/ for(i=j+1;i<=9;i++)for(j1=j+1;j1<=10;j1++)a[i][j1]-=(a[i][j]/a[j][j])*a[j][j1];} /*逐行消元，将矩阵转化为上三角矩阵*/ printf("x9为：%lf\n",x[9]=a[9][10]/a[9][9]);for(m=8;m>=1;m--){sum=0;for(n=m+1;n<=9;n++) sum+=x[n]*a[m][n]; printf("x%d 为：%lf\n",m,x[m]=(a[m][10]-sum)/a[m][m]);} /*从x9到x1逐个解出结果*/ }3.计算结果4.问题讨论（误差分析、上机出现情况等）（1）这道题目出现的最大问题在 for(j1=j+1;j1<=10;j1++)a[i][j1]-=(a[i][j]/a[j][j])*a[j][j1]; ，刚开始我错写为：for(j1=j;j1<=10;j1++)a[i][j1]-=(a[i][j]/a[j][j])*a[j][j1]; ，也就是j1应该从j+1开始取值，如果从j 开始取，会导致a[i][j1]-=(a[i][j]/a[j][j])*a[j][j1]; 一句中a[i][j]是零，这样后面的a[i][j1]全部都将是原来的数值，不会变化，导致错误。