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?。

?0289?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?

?0289?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?

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