2008级数值计算方法选题

（1）用Matlab语言或你熟悉的其他算法语言编程序。

（2）上机前充分准备，复习有关算法，写出计算步骤，反复检查程序。

（3）完成计算后写出实验报告，内容包括：所用算法语言，算法步骤叙述，变量说明，程序清单，输出计算结果，结果分析和小结等。

完成后请把所得结果写成论文形式完成，以PPT汇报，汇报结束修改后发到老师邮箱。

选题：

1． 编制以离散点{xi}i0的正交多项式{Pk(x)}为基的最小二乘拟合程序，并用于对下列数

据做三次多项式最小二乘拟合。

mxi yi 实验要求：

－1.0 －4.447 －0.5 －0.452 0.0 0.551 0.5 0.048 1.0 －0.447 1.5 0.549 2.0 4.552 取权(xi)1，求出拟合曲线yS(x)P(x)kxk，输出

kkk0k033k,Pk(x),kk(0及平方误差,||||22，并画出yS(x)的图形。

2． 给出积分 （1）I20xe2x2（2）dx，

234cotxdx，（3）I321dx。 2x1(0)(0)实验要求：（1）用龙贝格求积计算上述积分I的值。要求|Tk1Tk|106时结束，输出T表及I的近似值。（2）用5点高斯求积公式及复化3点高斯求积公式计算上述积分，并输出

（3）分析比较各种计算结果。 I的近似值。

3． 比较求一阶导数的数值方法。给定函数f(x)e,x[0.5,2]。

实验要求：利用某距离点函数值，必要时给定端点导数值，分别用中心差分，数值积分求导，三次样条求导和理查森外推计算f(x)的一阶导数，分析，比较各种方法的效果，说明精度与步长h的关系。

4． 给定方程组

1x3.016.031.99x11x14.161.23（1）1.272； 0.9874.819.34x31701x181032.09999962x5.9000012。 （2）5151x35x210214实验要求：

（1）用LU分解和列主元高斯消去法求解上述两个方程组。输出Ax=b中矩阵A及向量b，A=LU分解的L和U，detA及解向量x。

（2）将方程组（1）中系数3.01改为3.00，0.987改为0.990。用列主元高斯消去法求解，输出列主元行交换次序、解向量x及detA，并与（1）中结果比较。

（3）将方程组（2）中系数2.099999改为2.1，5.900001改为5.9。用列主元高斯消去法求解，输出解向量x及detA，并与（1）中结果比较。 （4）对结果做分析。

5． 用迭代方法求解线性方程组。用迭代法求解方程组Ax = b，其中A∈ R1010为五对角

矩阵，它的每条对角线元素是常数，设为

31214A实验要求：

12312141412312141412312141412312 14123（1）选取不同的初始向量x(0)和不同的方程组右端项向量b，给定迭代误差要求，分别用Jacobi、Gauss-Seidel以及SOR（松弛系数取12的不同值） 迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛，记录迭代次数。比较并分析实验结果， 要求迭代误差满足x(k1)x(k)105。

（2）取定右端向量b 和初始向量x(0)（自己定），将A的主对角线元素成倍增长若干次（次数自己定），非主对角线元素不变，重复（1）中过程，比较收敛速度和实验结果，迭代误

差同（1）。

6． 求非线性方程及方程组的根，准确到106。给定方程分别为：

23x12x20(())T x(1,1)（1）x3x2e0； （2）。 233x1x2x1102x实验要求：

（1）用你自己设计的一种线性收敛迭代法求方程（1）的根，然后再用斯蒂芬森加速迭代计算。

（2）用牛顿法求方程（1）的根，输出迭代初值，各次迭代值及迭代次数，并与（1）的结果比较。

（3）用牛顿法求（2）的根，输出迭代次数及解向量x的近似。

7． 用QR算法求矩阵特征值：

24621（1）A231； （2）Η011100345645673678。

020010实验要求：

（1）根据QR算法原理编制求（1）及（2）中矩阵全部特征值的程序并输出计算结果（要求误差105）。

（2）直接用现有数学软件求（1），（2）的全部特征值，并与（1）的结果比较。

8． 求初值问题的数值解，给定初值问题为：

y50y50x22x,0x1,1yy2,1x2,（1） （2） xx1y(0).y(1)1.3实验要求：

（1）用改进的欧拉法（取h0.05）四阶R-K方法（取h0.1）求（1）的数值解，并输出xi10.1i(i0,1,,10)的数值解yi。

（2）用经典四阶R-K方法解（2），步长h分别取为h0.1， 0.025， 0.01计算，并输出

xi0.1i(i0,1,,10)各点的数值解yi及准确解y(xi)，并分析结果。（初值问题（2）的

准确解y(x)150xex2） 3

9．常微分方程初值问题

y'ycos2x2sin2x2xex 0（1）分别取步长h = 0.1，0.01，0.001，采用改进的Euler方法、4阶经典R－K方法和4阶Adams预测-校正方法计算初值问题。以表格形式列出10个等距节点上的计算值

和精确值，并比较他们的计算精确度。其中多步法需要的初值由经典R-K法提供。运算时，取足以表示计算精度的有效位。

(2) 取h0.001，仍然用这两种方法（4阶经典R－K方法和4阶Adams预测-校正方法）计算(不必保留中间结果.只需显示最后一步的计算值)，比较两种方法所花的计算机CPU时间。

2xcos(2x)。

2'2sin(lnx),1x2y\"y2y10．（1）  xxx2y11,y22其精确解为y(x)c1xc231sin(lnx)cos(lnx),其中 210x101(812sin(lin2)4cos(ln2))0.039207013270 11c1c21.1392070310c2y\"xy'3y4x,0x13（2），其精确解为yxxx

y(0)0,y(1)2实验要求：对每个边值问题，分别取h=0.2和h=0.1,用差分法计算,并利用两个步长计算得到的近似解进行外推。列出计算解，外推解和精确解，并进行比较。

11、线性方程组的稳定性。体会用直接方法和迭代方法求解病态线性方程组时，各种方法的收敛性问题。

考虑方程组Hx = b 的求解，其中系数矩阵

H 为Hilbert病态矩阵：

H(hij)nn,hij试验要求：

1,ij1i,j1,2,,n

（1）取定右端向量b（自己定，可通过先给定解如各分量均为1，再计算出右端b），选择

方程组的阶数为n = 5，分别用LU 分解方法、Jacobi、Gauss-Seidel 和 SOR 迭代方法求解方程组，各自的结果如何？将计算结果与问题的解比较，结论如何？

(2)逐步增大问题的维数（n = 10、20、30…），计算得结果如何？计算结果说明了什么？

12、 矩阵特征值和特征向量的数值计算方法 考虑Hilbert 矩阵

H(hij)nn,hij试验要求：

1,ij1i,j1,2,,n

（1）取n = 5，分别编写实用的QR迭代算法（不是用Matlab 中提供的函数“eig”）、隐式对称QR 迭代算法、Jacobi 方法计算该矩阵的全部特征值，并用带位移的幂法求所有特征向量，比较并分析实验结果。

（2）分别利用原理型的QR算法和Matlab中提供的函数“eig”计算下面矩阵的特征值并比较其结果。

24A000

345645673678

02001013、高维积分数值计算的

Monte-Carlo方法

相关原理：设D为s维空间Rs的一个区域，f(x)DRsR，考虑高维积分

If(x)dx

D高维积分一般采用随机抽样方法，即Monte-Carlo方法．

Monte-Carlo方法有以下两种形式：

１．在D上随机抽取n个点x(1),x(2),,x(n)，计算f(x)在这些点处的函数值

f(x(1)),f(x(2)),,f(x(n))

计算f(x)在D上的样本均值

1nˆfnf(x(i))

ni0ˆ自然，我们希望fnI，m(D)为D的测度（面积、体积等），因而 m(D)m(D)nf(x)dxf(x(i)) ni0If(x,y)dxdy

D ID２．方法２以计算二重积分为例说明其基本思想，设

我们总可以将积分域D包含在矩形区域[a,b][c,d]中，而且不失一般性假设．因此，被积函数下面的柱体0f(x,y)M中（f(x,y)0时的处理方式相同）

可以包含在立方体W[a,b][c,d][0,M]中．在W中随机选取一个点其中x,y,z分别是在区间[a,b],[c,d],[0,M]中均匀分布的随机数．根K(x,y,z)，

据几何概率的理论，点K落在柱体中的概率等于的体积与W的体积之比，

即

P(KKW)所以

IM(ba)(dc)P(KKW)

计算时用抽样的方法算出K落在中的频率代替概率，得到

If(x,y)dxdyM(ba)(dc)DI

M(ba)(dc)n N其中N是抽样总数，n是抽样点K落在之中的频数．

对于一般的区域D，计算其测度（只要理解为平面上的面积和空间中的体积）的一般方法是：先找一个规则的区域A包含D，且A的测度m(A)是已知的．生成A中N个均匀分布的点x(1),x(2),,x(N)，如果其中有n个落在区域D中，则区域D的测度近似为

m(D)nm(A) N函数f(x)在区域D上的积分可以近似为区域D的测度与函数f(x)在区域D中的n个随机值的平均值的乘积，即

Df(x)dxm(D)1f(x(k)) nx(k)Dm(A)f(x(k)) Nx(k)D实验要求

１．假设冰淇淋的下部为一锥体而上面为一半球，考虑冰淇淋锥的体积问题：计算锥面z2x2y2上方和球面x2y2(z1)21的内部区域的体积．如果使用球面坐标，该区域可以表示为如下的积分：

2cos204002sinddd

用Monte-Carlo方法计算该积分；

２．显然这样的冰淇淋可以装在如下立方体的盒子里 1x1,1y1,0z2 而该立方体的体积为８．生成这个盒子里均匀分布的随机点，落入冰淇淋锥中的点数与总点数之比再乘以８得到的就是冰淇淋锥的体积．将该结果与１中的结果进行比较；

３．类似的方法可以计算复杂区域的测度（面积或体积）．试求由下列关系所界定区域的测度：

0x1,1y2,1z3exy（１）

ysinz01x3,1y4（２）x3y329

yex20x1,0y1,0z1x2sinyz（３）

xzey1