《计算方法》复习题-1.doc

第一章 绪 论 43 第一章 绪 论复习题例1计算，取，采用下列算式计算：（1）；（2）；（3）；（4）。问哪一个得到的结果最好？解 显然所以，这4个算式是恒等的，但当取计算时，因为（2），（3）都涉及到两个相近数相减，使有效数字丢失，而（1）在分母算式上的乘幂数比算式（4）大，所以算式（4）最好。事实上，当取时，有，再由的误差也可直接估计出每个算式的误差，显然，算式（4）误差最小。具体计算可得：（1）；（2）；（3）；（4）。比较可得用第（4）个算式所得的结果更接近于。例1.8 建立积分的递推关系式，并研究它的误差传递。解 由和可建立下列递推公式 (*)计算出后，由递推关系式可逐次求出的值。但在计算时有舍入误差，因此在使用递推关系式中，实际算出的都是近似值。即现在来研究误差是如何传递的。设有误差，假设计算过程中不产生新的舍入误差，则由(*)式可得从而有即原始数据的误差对第步的影响使该误差扩大了倍。当较大时，误差将淹没真值，因此递推公式(*)是数值不稳定的。现在从另一方向使用这一公式 (*)只要给出的一个近似值，即可递推得到，类似于上面的推导可得每递推一步误差缩小到原值的，所以递推公式(*)是数值稳定的。由于时，所以有估计式于是取可得另一算法：由此可见，对于同一数学问题，使用的算法不同，效率也大不相同，只有选用数值稳定性好的算法，才能求得较准确的结果。基于Mathematica的数值计算实例例1 计算有位有效数字的近似值，并列表。解 Mathematica程序：TableNE,n,NPi,n,n,1,10;TableForm%运行结果： 3. 3.2.7 3.12.72 3.142.718 3.1422.7183 3.14162.71828 3.141592.718282 3.1415932.7182818 3.14159272.71828183 3.141592652.718281828 3.141592654例2 用程序计算有位有效数字的近似值。解 Mathematica程序：TableNSqrt1500,n,N12(1/6),n,n,10,15;TableForm%运行结果：38.72983346 1.51308574938.729833462 1.513085749438.7298334621 1.5130857494238.72983346207 1.51308574942338.729833462074 1.513085749422938.7298334620742 1.5130857494229例3 计算的近似值。解 Mathematica程序：Cos75.5 Degree,ArcTan34.7 Degree,Log5,79/N;TableForm%运行结果：0.250380.5445482.71489例4 二项式系数定义为，利用该定义计算。解 Mathematica程序：CCn_,k_:=n!/(k!*(n-k)!)CC50,36运行结果：937845656300例5 分别给出前20个素数及第100个素数。解 Mathematica程序：TablePrimen, n, 1, 20Prime100运行结果：前20个素数：2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71第100个素数541例6 用秦九韶法计算在处的值，并验证之。解 Mathematica程序：ak_:=k+1;s8=a8;sk_:=x*sk+1+ak;x=2;s0运行结果：4097。第二章 插值与拟合复习题例2.1 插值函数作为被插函数的逼近，可以用作函数值的近似计算。已知，构造二次拉格朗日插值多项式。（1）计算；（2）估计误差并与实际误差相比较。解 （1）以插值点(27,3), (64,4), (125,5)代入插值公式，得=(2) 由误差公式有记在27，125上是单调递减函数。实际误差：。 例2.2已知sin0.320.314 567，sin0.340.333 487，sin0.360.352 274，用线性插值及抛物插值计算sin0.336 7的值并估计截断误差。分析题目中相当于告诉了插值条件。考虑到0.336 7位于0.32与0.34之间，根据插值法的特点，线性插值时，取0.32和0.34作为插值节点；抛物插值时，三个点全取。由于一次、二次插值函数表达式较简单，可采用牛顿型公式，误差估计用拉格朗日型余项表达式。解用线性插值计算，则sin0.336 70.314 567其截断误差限为其中，于是0.92用抛物插值计算，有0.0167其截断误差限为其中，于是0.02330.203注=0.330374与具有六位有效数字的正弦函数表完全一样，说明用二次插值精度已相当高了。例2.3 用插值点（2，4），（3，9），（5，25）分别构造拉格朗日插值函数和牛顿插值函数，并计算和。解（1）以插值点（2，4），（3，9），（5，25）代入插值公式，得+=代入可得。（2）做出插值点（2, 4）（3, 9）（5, 25）的差商表：024139(9-4)/(3-2)=52525(25-9)/(5-2)=8(8-5)/(5-2)=1代入可得。例2.4 设，和取值如下：-1-0.500.514.22.451.20.450.2分别构造，并比较结果。解 以（-1, 4.2），（0, 1.2）和（1, 0.2）为插值点，则 以（-1, 4.2），（-0.5, 0.45），（0, 1.2）和（1, 0.2）为插值点，则以（-1, 4.2），（-0.5, 2.45），（0, 1.2），（0.5, 0.45）和（1, 0.2）为插值点，则注意到，所以时，从而有。例2.5已知函数的函数表如下：0.400.550.650.800.901.050.410 750.578 150.696 750.888 111.026 521.253 82求四次牛顿插值多项式，并由此求的近似值。分析表中给出六对数据，故最高可构造五次多项式。但由于0.596接近于，因此可取前五对数据来做差商表。解构造差商表如下：一阶差商二阶差商三阶差商四阶差商0.400.550.650.800.900.410 750.578 150.696 750.888 111.026 521.116 001.186 001.275 731.384 100.280 000.358 930.433 480.197 330.213 000.031 34故四次牛顿插值多项式为 于是0.631 95。例2.6利用差分性质证明：123n分析 本题结论是知道的。但这里要用差分性质来证，那就必须对某个函数构造差分。因此关键要寻找一个函数，使其利用差分的某个性质，经过推导，得出题断。下面仅对的结论给出两种证法，的结论类似地可证。证法一容易看出，当令12n时，为此可得差分表如下：012011+21211100利用差分公式，取0，则 即123n证法二由于 +所以定义函数，只要证明（1，2，n1）且即可证得结论。事实上，（1，2，n1），是显然的。故123n例2.7 试对下列数据做出形如的拟合曲线。2357816224661解 写出的矛盾方程组：法方程为：解得：，。所以有例2.8 试给出下列数据0.30.50.60.70.91.377311.487661.538791.586531.67的形如的拟合函数。解 矛盾方程为：法方程组为：解方程组得：，。所以有例2.9 给出下列数据1.251.371.451.691.7720.866624.481928.66739.672646.22试对数据做出形如的拟合函数。解 对函数两边取自然对数：得到数据表：1.251.371.451.691.773.038153.197933.355753.68.663.83341得矛盾方程组：法方程组为：解方程组可得：，。所以有基于Mathematica的数值计算实例例1 已知求的二次插值多项式。解 Mathematica程序：A1 = 0, 0, 1, 3, 2, 9;f = InterpolatingPolynomialA1, x;Expand%运行结果： 例2 已知函数表12345141310711求出Lagrange插值多项式，并计算处的的近似值。解 Mathematica程序：r0x_,x0_,a_,b_,c_,d_:=(x-a)(x-b)(x-c)(x-d)/(x0-a)(x0-b)(x0-c)(x0-d)L4x_:=r0x,1,2,3,4,5*14+r0x,2,1,3,4,5*13+r0x,3,1,2,4,5*10+r0x,4,1,2,3,5*7+r0x,5,1,2,3,4*11L4x/NSimplify%L31.5/N运行结果：插值多项式为0.583333 (-5. + x) (-4. + x) (-3. + x) (-2. + x) - 2.16667 (-5. + x) (-4. + x) (-3. + x) (-1. + x) + 2.5 (-5. + x) (-4. + x) (-2. + x) (-1. + x) - 1.16667 (-5. + x) (-3. + x) (-2. + x) (-1. + x) +0.458333 (-4. + x) (-3. + x) (-2. + x) (-1. + x)插值多项式的简式为插值为。例3 已知14个点的坐标如下表，根据表中数据构造一插值多项式，并由此计算点的函数值的近似值（不要求给出插值多项式）。解 已知数据表：11.52.02.53.03.54.04.55.02.718284.481697.3890612.182520.085533.115554.598254.598290.0171Mathematica程序：A=Tablex,Expx,x,1,5,0.5/Ng1=ListPlotTableA,Prolog-AbsolutePointSize6, DisplayFunction-Identity;InterpolationA,InterpolationOrder-1g2=Plot%x,x,1,5,PlotStyle-Thickness0.005,DisplayFunction-IdentityShowg1,g2,DisplayFunction-$DisplayFunction%4.25运行结果：所求插值。例4 已知函数的函数表如下00.10.20.30.40.50.600.099830.198670.295520.389420.4794260.56464(1) 分别用线性插值、二次插值和三次插值求的近似值，并估计截断误差；(2) 试用4次等距节点插值公式计算的近似值，并估计误差。解 Mathematica语句：Tablea = 0.3,0.4, 0.5, 0.6, 0.29552, 0.38942, 0.47943, 0.56464;Transposea2, 3;y1 = InterpolatingPolynomial%, x / ExpandTransposea2, 3, 4;y2 = InterpolatingPolynomial%, x / ExpandTransposea1, 2, 3, 4;y3 = InterpolatingPolynomial%, x / Expand Thickness0.01, Thickness0.008, Thickness0.005, PlotLegend - y1, y2, y3, LegendPosition - 1, 0MapNumberForm#, 10 &, Sin0. 48, y1, y2, y3 /. x - 0.48运行结果： 一次插值多项式 二次插值多项式 三次插值多项式下面是的真值与三种插值的比较：0.4617791755, 0.461428, 0.461812, 0.46178288 最后估计误差的Mathematica程序为：resErrorx_, n_, xp_List, k_ := DSint, t, n + 1/(n + 1)!)* Product(x - xp)j + 1 /. t - k, j, 0, nresError0.48, 1, 0.4, 0.5, 0.5 /. t - 0.5resError0.48, 2, 0.4, 0.5, 0.6, 0.4 /. t - 0.4resError0.48, 3, 0.3, 0.4, 0.5, 0.6, 0.3 /. t - 0.3以下是三种插值在计算时的误差的运行结果0.00038354-0.000029474(2) Mathematica程序：xk=k*0.1;y0=0;y1=0.09983;y2=0.19867;y3=0.29552;y4=0.38942;y5=0.479426;y6=0.56464;f1i_:=yi+1-yif2i_:=f1i+1-f1if3i_:=f2i+1-f2if4i_:=f3i+1-f3if5i_:=f4i+1-f4if6i_:=f5i+1-f5iA=y0,y1,y2,y3,y4,y5,y6,0,f10,f11,f12,f13,f14,f15, 0,0,f20,f21,f22,f23,f24,0,0,0,f30,f31,f32,f33, 0,0,0,0,f40,f41,f42,0,0,0,0,0,f50,f51,0,0,0,0,0,0,f60;TransposeA/N;NMatrixForm%,4a0=y0;a1=f10;a2=f20;a3=f30;a4=f40;a5=f50;a6=f60;NNx=Sumak*Product(t-j),j,0,k-1,k,0,6/NNExpand%,3/Chop%/.t-0.48运行结果：所得差分表为：0 0 0 0 0 0 00.09983 0.09983 0 0 0 0 00.1987 0.09884 -0.00099 0 0 0 00.2955 0.09685 -0.00199 -0.001 0 0 00.3894 0.0939 -0.00295 -0.00096 0.00004 0 00.4794 0.09001 -0.003894 -0.000944 0.000016 -0.000024 00.5646 0.08521 -0.004792 -0.000898 0.000046 0.00003 0.000054插值多项式为：0.09983 t - 0.00099 (-1. + t) t - 0.001 (-2. + t) (-1. + t) t + 0.00004 (-3. + t) (-2. + t) (-1. + t) t -0.000024 (-4. + t) (-3. + t) (-2. + t) (-1. + t) t + 0.000054 (-5. + t) (-4. + t) (-3. + t) (-2. + t) (-1. + t) t化简得： 则0.468457。例5 过两点构造一个三次插值多项式，满足条件。解 Mathematica程序：x0=0;x1=1;y0=0;y1=1;y0=2;y1=2;wx_,1_:=Product(x-xk),k,0,1;lx_,k_:=wx,1/(x-xk)(Dwx,1,x/.x-xk)Dlx_,k_:=Dlx,k,x/.x-xkhx_,k_:=(1-2(x-xk)(Dlx,k/.x-xk)(lx,k2);Hx_,k_:=(x-xk)(lx,k2);Tablehx,k,k,0,1;TableHx,k,k,0,1;HHx_,n_:=Sumyk*hx,k+ykHx,k,k,0,n;HHx,1Expand%运行结果： 化简为： 例6 求数据表12345678910-1-0.8-0.6-0.40.200.20.40.60.8-0.32090.42950.98261.23922.00002.44453.03343.56013.77874.1564的最小二乘二次拟合多项式。解 最小二乘求解的Mathematica程序：ClearxX=Table-1+0.2*i,i,0,9;Y=-0.3209,0.4295,0.9826,1.2392,2.0000,2.4445,3.0334,3.5601,3.7787,4.1564;qa_,b_,c_:=Sum(a+b*Xk+c*Xk2-Yk)2,k,1,10SolveDqa,b,c,a=0,Dqa,b,c,b=0,Dqa,b,c,c=0,a,b,c运行结果：拟合多项式系数a - 2.49261, b - 2.42838, c - -0.351241拟合多项式2.49261 + 2.42838 + 0.351241LL=TableXk,Yk,k,1,9g1=ListPlotLL,Prolog-AbsolutePointSize10f=FitLL,1,x,x2,xg2=Plotf,x,-1,1Showg1,g2运行得拟合效果：例7（双曲拟合）给定数据x1.01.11.21.31.41.5y0.9710.7720.5970.4290.1680.065求形如的拟合函数。解 Mathematica程序：ClearX,Y,f,LL,g1,g2L=1.0,0.971,1.1,0.772,1.2,0.597,1.3,0.429,1.4,0.168,1.5,0.065;X=1.0,1.1,1.2,1.3,1.4,1.5;Y=0.971,0.772,0.597,0.429,0.168,0.065;YY=1/Y;LL=TableXn,YYn,n,1,6f=FitLL,1,x,xg1=ListPlotL,Prolog-AbsolutePointSize15g2=Plot1/f,x,0,3Showg1,g2运行得到结果： ，。从而得到拟合曲线方程为 拟合曲线为 例9 设有一形为的函数，现测得函数值表如下：0.10.20.30.40.50.60.70.80.90.340.570.891.341.491.802.553.485.00试用最小二乘法确定和。解 Mathematica程序：ClearA,X,Y,LL,f,g1,g2X=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9; Y=0.34,0.57,0.89,1.34,1.49,1.80,2.55,3.48,5.00;L=TableXn,Yn,n,1,9LL=TableXn,LogYn,n,1,9f=FitLL,1,x,x运行结果为：-1.18014 + 3.0968 x。即A=-1.18014;b=3.0968;a=ExpA;y=a*Expb*xg1=ListPlotL,Prolog-AbsolutePointSize5g2=Ploty,x,0,0.9Showg1,g2运行得结果为：则原函数可近似为： 拟合曲线为：第三章 线性方程组的解法复习题例3.1 用高斯消元法解方程组解 取小数位4位对做高斯消元：=由同解方程组作回代过程有： 例3.2 用列主元素法解方程组：解 对做选主元及消元过程：由同解方程组作回代过程有：例3.3 分别用Doolittle和Crout分解法求解方程组解 方程系数矩阵，常数项。(1) Doolittle分解：对等式按的第一行、的第一列，按的第二行、的第二列，的第三行的次序，根据矩阵乘法与矩阵相等，对比等号两边矩阵元素，易得矩阵,元素：，求解，得。求解，得。(2) Crout分解对等式 根据矩阵乘法与矩阵相等，对比等号两边矩阵元素，易得矩阵元素：，求解，得。求解，得。例3.4用矩阵的直接三角分解法解方程组分析这是常规计算题，只需按矩阵的三角分解过程认真计算即可。解设由矩阵乘法可逐行，逐列分别求出和：有5，3，6，4再解上三角方程组得2，2，1，1。注由于三角分解过程十分规律，也可按紧凑格式直接得到其中折线右上部的元素（竖线左侧）为上三角阵的元素，最后一列为的元素，故，例3.5用平方根法（Cholesky分解）解方程组=分析由于系数矩阵对称正定，故一定有分解形式，其中为下三角阵，然后由矩阵乘法即可求出的元素。解设由矩阵乘法得，由解得，。再由得，。例3.11 计算向量的1范数、2范数和无穷范数。解 由向量范数定义得 =1+2+1.5=4.5 =2.69258 =2例3.12 计算矩阵的1范数、2范数、范数和Frobenius范数。解 由矩阵范数定义： 1=max1.1+2.5，+3.5=5.5 =max1.1+，2.5+3.5=6 F 特征值123.6541，20.05591 用Mathematica编写的程序： A=1.1，-2，2.5，-3.5;MatrixForm%;A1=MaxSumAbsAn，n，1，2B=TransposeA;AI=MaxSumAbsBm，m，1，2AF=SqrtSumAi，j2，i，1，2，j，1，2/NT=TransposeA.A ;MatrixForm%; EigenvaluesA; A2=NSqrtMaxEigenvaluesT，10运行结果为：5.56.4.869294.86354706例3.16 用雅可比方法求解方程组方程的准确解为。解 (1) 由系数矩阵构造雅可比迭代矩阵：，取。计算结果见下表： 1234567891.5 0.5 -8 7.59.83333 -9 -4.75 9.515.3056 -39.0833 10.6667 30.083313.5648 -45.5556 80.4028 69.7361-55.1744 31.1435 98.3241 76.6991-144.488 324.022 -24.6844 292.878-178.648 558.268 -667.566 642.881363.378 52.0271 -1219.02 551.4571444.93 -2667.53 -495.445 2719.56可以看到迭代发散。的特征值。所以有(2) 将方程的次序调换可以得到构造迭代格式：写成迭代矩阵形式：迭代收敛，取。结果见下表：k 1234567891011121.75 -1.6 1 0.751.9 -1.9 0.95 0.31.9625 -1.94 1 0.06251.985 -1.985 0.9925 0.0451.99438 -1.991 1 0.0093751.99775 -1.99775 0.998875 1.99916 -1.99865 1 1.99966 -1.99966 0.999831 1.99987 -1.9998 1 1.99995 -1.99995 0.999975 1.99998 -1.99997 1 1.99999 -1.99999 0.999996 例3.17 用高斯塞德尔方法求解方程组其中方程的准确解为。解 高斯塞德尔迭代格式： 写成迭代矩阵形式： 因为。所以，高斯塞德尔迭代格式收敛速度比雅可比迭代快。 取为初始值，下表列出迭代计算结果：k123456781.251.96252.018132.004572.000231.999891.999972.-1.7-2.045-2.01825-2.00185-1.99973-1.99989-1.99999-21.151.02751.000040.9990910.9998320.9999991.0000112.70.71250.0556250.01640420.004338723.4369510-49.96610-52.6418910-6Mathematica程序：A=4，1，-1，2，5，2，1，1，3;MatixForm%;b=5，-4，3;I1=IdentityMatix3;DD=A1，1，0，0，0，A2，2，0，0，0，A3，3;L=-0，0，0，A2，1，0，0，A3，1，A3，2，0;DDLN=InverseDD-L;G=I1-DDLN.A;MatrixForm%;p=MaxAbsEigenvaluesNGR=p-1f=DDLN.b;x0=1，1，1;xn_：=G.xn-1+f;Tablexn，n，1，8;MatrixForm%/NTableMaxAbsxn-xn-1，n，1，8;MatrixForm%/NLinearSolveA，b例3.18 如何对方程组进行调整，使得用高斯-塞德尔方法求解时收敛？并取初始向量，用该方法求近似解，使。 解 将第三个方程调到第一行后有这是主对角线严格占优方程组，故用高斯-塞德尔迭代法求解一定收敛。迭代格式为由，得因为 ，故取最后结果为例3.19 给定方程组 证明：雅可比方法发散，而高斯-塞德尔收敛。 证明 对雅可比方法，迭代矩阵为设其特征值为，则，故雅可比方法发散。对高斯-塞德尔方法，迭代矩阵为=显然，其特征值为，故高斯-塞德尔收敛。例3.20 设为正交矩阵，求证：线性方程组用高斯-塞德尔方法求解必收敛。证明：因为正交矩阵，故，设是的特征值，则有，使两边与作内积，有从而有，又的特征值为，故的特征值不为零，从而非奇异。故正定，所以高斯-塞德尔方法收敛。例3.21 对方程组(1) 给出解方程组的雅可比迭代矩阵，并讨论迭代收敛条件；(2) 给出解方程组的高斯-塞德尔迭代矩阵，并讨论迭代收敛条件。解 (1) 雅可比迭代矩阵： 其特征多项式为：故的谱半径：即雅可比迭代收敛的充要条件为。(2) 高斯-塞德尔迭代矩阵：的特征多项式为：故的谱半径：即高斯-塞德尔迭代收敛的充要条件为。基于Mathematica的数值计算实例例1 求解线性方程组 解 Mathematica程序： DetA0， A=3,4,-5,7,2,-3,3,-2,4,11,-13,16,7,-2,1,3;MatrixForm%DetA r=SumRowReduceAi,i,i,1,4 NullSpaceA运行结果为：行列式等于0，系数矩阵的秩为2，有非零解，基础解系为-13, -20, 0, 17, 3, 19, 17, 0。可得方程组的通解为：例2 求解线性方程组 解 Mathematica程序：A=-1,3,2,-1,3,-5,-3,2,1,3,7,5,7,-3,-5,-1;MatrixForm%DetA NullSpaceA运行结果为：行列式的值为，方程组只有零解。例3 求解线性方程组 解 Mathematica程序：ClearA,bA=2,1,-4,-1,6,-2,-2,2,0,5,14,3,4,-3,-10,-1;MatrixForm%;b=8,4,-4,8;DetA r=SumRowReduceAi,i,i,1,4 NullSpaceA LinearSolveA,b运行结果为：行列式等于0，系数矩阵的秩为3，有非零解，基础解系为-1, 1, -1, 3，一个特解为1, 2, -1, 0。则原方程组的通解为：例4 方阵由前36个相邻整数构成，求的特征多项式、特征方程、特征值与特征向量。解 Mathematica程序：A=PartitionRange36, 6;MatrixForm%CharacteristicPolynomialA,CharacteristicPolynomialA,0EigenvaluesNA/ChopEigenvectorsA/N;MatrixForm%运行结果： 特征多项式 特征方程116.412, -5.41182, 0, 0, 0, 0 特征根 特征向量例5 用高斯-若当(Gauss-Jordan)消元法解线性方程组解 Mathematica程序：A=2,3,1,6,1,4,-1,4,1,-2,1,1;MatrixForm%RowReduceA;MatrixForm%运行结果：则求得原方程组的解为：，。例6 设向量，求向量范数。解 Mathematica程序： x=2,-5,1,3,0;=SumAbsxi,i,1,5 运行结果：11=Sqrtx.x 运行结果： =MaxTableAbsxi,i,1,5 运行结果：5例7 设矩阵求矩阵的范数。解 Mathematica程序： A=3,2,0,1,-1,0,-1,4,2,1,1,3,2,3,1,5;MatrixForm%;A1=MaxSumAbsAn,n,1,4A100=MaxTableSumAbsAn,m,m,1,4,n,1,4T=TransposeA.A;MatrixForm%;EigenvaluesNT;A2=NSqrtMaxEigenvaluesNT 运行结果： 13 11 8.25543例8 考察阶希尔伯特(Hilbert)矩阵的性态解 Mathematica通用程序：ClearBB=Table1/(i+j),i,1,k,j,0,k-1;MatrixForm%BN=InverseB;B1=MaxSumAbsBn,n,1,k;BN1=MaxSumAbsBNn,n,1,k;condA1=B1*BN1/NB100=MaxTableSumAbsBn,m,m,1,k,n,1,k;BN100=MaxTableSumAbsBNn,m,m,1,k,n,1,k;condA100=B100*BN100/N当取4时, 运行结果： 28375 28375例9 用Jacobi迭代法求解方程组 解 Mathematica程序：A = 7, -4, 2, 4, 10, -1, 6, 3, 9;MatrixForm%;b = 18, 40, 29;I1 = IdentityMatrix3;DD = A1, 1, 0, 0, 0, A2, 2, 0, 0, 0, A3, 3DDN = InverseDD;J = I1 - DDN.A;MatrixForm