实验四 数值积分与数值微分

实验4.1 （高斯数值积分方法用于积分方程求解）

问题提出：线性的积分方程的数值求解，可以被转化为线性代数方程组的求解问题。而线性代数方程组所含未知数的个数，与用来离散积分的数值方法的节点个数相同。在节点数相同的前提下，高斯数值积分方法有较高的代数精度，用它通常会得到较好的结果。

实验内容：求解第二类Fredholm积分方程

y(t)??k(t,s)y(s)ds?f(t),a?t?b

ab首先将积分区间[a,b]等分成n份，在每个子区间上离散方程中的积分就得到线性代数方程组。

实验要求：分别使用如下方法，离散积分方程中的积分 1.复化梯形方法；2.复化辛甫森方法；3.复化高斯方法。求解如下的积分方程。

21tttey(s)ds?ee1）y(t)?，方程的准确解为； ?0e?1112）y(t)??esty(s)ds?1?e2t?et，方程的准确解为et；

02 比较各算法的计算量和误差以分析它们的优劣。

实验4.2（高维积分数值计算的蒙特卡罗方法）

问题提出：高维空间中的积分，如果维数不很高且积分区域是规则的或者能等价地写成多重积分的形式，可以用一元函数积分的数值方法来计算高维空间的积分。蒙特卡罗方法对计算复杂区域甚至不连通的区域上的积分并没有特殊的困难。

实验内容：对于一般的区域?，计算其测度（只要理解为平面上的面积或空间中的体积）的一般方法是：先找一个规则的区域A包含?，且A的测度是已知的。生成区域A中m个均匀分布的随机点pi,i?1,2,?,m，如果其中有n个落在区域?中，则区域?的测度m(?)为n/m。函数f(x)在区域?上的积分可以近似为：区域?的测度与函数f(x)在?中n个随机点上平均值的乘积，即

??f(x)dx?m(?)?1f(pk) ?npk?? 实验要求：假设冰琪淋的下部为一锥体而上面为一半球，考虑冰琪淋体积问

题：计算锥面z2?x2?y2上方和球面x2?y2?(z?1)2?1内部区域的体积。如果使用球面坐标，该区域可以表示为如下的积分：

?/42cos?02???0?0?2sin?d?d?d?

用蒙特卡罗方法可以计算该积分。

另一方面，显然这样的冰琪淋可以装在如下立方体的盒子里

?1?x?1,?1?y?1,0?z?2

9

而该立方体的体积为8。只要生成这个盒子里均匀分布的随机点，落入冰琪淋锥点的个数与总点数之比再乘以8就是冰琪淋锥的体积。比较两种方法所得到的结果。

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

?0?x?1,1?y?2,?1?z?3?(1)?ex?y ?sin(z)y?0??1?x?3,?1?y?4?(2)?x3?y3?29 ?y?ex?2??0?x?1,0?y?1,0?z?1?(3)?x2?siny?z ?x?z?ex?1?

相关MATLAB函数提示： diff(x) 如果x是向量，返回向量x的差分；如果x是矩阵，则按各列作差分 diff(x,k) k阶差分 q=polyder(p) 求得由向量p表示的多项式导函数的向量表示q Fx=gradient(F,x) 返回向量F表示的一元函数沿x方向的导函数F'(x)，其中x是与F同维数的向量 z=trapz(x,y) x表示积分区间的离散化向量；y是与x同维数的向量，表示被积函数；z返回积分的近似值 z=guad(fun,a,b,tol) 自适应步长Simpson积分法求得Fun在区间[a,b]上的定积分，Fun为M文件函数句柄，tol为积分精度 z=dblquad(fun,a,b,c,d,tol,method) 求得二元函数Fun(x,y)的重积分 z=triplequad(fun,a,b,c,d,e,f,tol,method) 求得三元函数Fun(x,y,z)的重积分 10

实验五 解线性方程组的直接方法

实验5.1 （主元的选取与算法的稳定性）

问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。

实验内容：考虑线性方程组

Ax?b,A?Rn?n,b?Rn

编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。

实验要求：

?61??7??861??15?????（1）取矩阵A??????,b????，则方程有解x*?(1,1,?,1)T。

????861???15???86????14??取n=10计算矩阵的条件数。让程序自动选取主元(顺序消元)，结果如何？

（2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。

（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。

（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。

实验5.2（线性代数方程组的性态与条件数的估计） 问题提出：理论上，线性代数方程组Ax?b的摄动满足

?xcon(dA)?x1?A?1?A??A?b????A?b? ??矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算A?1通常要比求解方程Ax?b还困难。

实验内容：MATLAB中提供有函数“condest”可以用来估计矩阵的条件数，

它给出的是按1-范数的条件数。首先构造非奇异矩阵A和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动?A和?b，使得?A和?b充

11

分小。

实验要求：

??b??b，以1-范数，（1）假设方程Ax=b的解为x，求解方程(A??A)x给出

?xx???xxx的计算结果。

（2）选择一系列维数递增的矩阵（可以是随机生成的），比较函数“condest”

所需机器时间的差别.考虑若干逆是已知的矩阵，借助函数“eig”很容易给出cond2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比较。

（3）利用“condest”给出矩阵A条件数的估计，针对（1）中的结果给出

?xx的理论估计，并将它与（1）给出的计算结果进行比较，分析所得结果。注意，如果给出了cond(A)和A的估计，马上就可以给出A?1的估计。 （4）估计著名的Hilbert矩阵的条件数。

H?(hi,j)n?n,hi,j?1,i?j?1i,j?1,2,?,n

思考题一：（Vadermonde矩阵）设

?1??1 A??1????1?x0x1x2?xn2x0?x12?2x2???2xn??ni???x0?0n?i???x0n?xi?n?1x1????i?0n?x2，b??ni?，

??x2?????i?0?n????xn??ni???xn??i?0?其中，xk?1?0.1k,k?0,1,?,n，

（1）对n=2,5,8，计算A的条件数；随n增大，矩阵性态如何变化？

（2）对n=5，解方程组Ax=b；设A的最后一个元素有扰动10-4，再求解Ax=b （3）计算（2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。 （4）你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗？

12