数值实验[综合]

数值实验一

数值实验综述：本章主要介绍数值分析的基本知识，其重点是对数值方法的评价与构造算法应注意的问题。本章数值实验的中心就是通过具体的数值实验来进一步熟悉MATLAB 的使用和理解算法的特点。如算法的计算复杂性，截断误差分析与收敛性，误差传播与算法的稳定性等。

实验1.1 病态问题

算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化比较敏感，反之属于好问题。

数值分析的大部分内容，如解线性代数方程组、求解矩阵特征值问题、解非线性方程及方程组等都存在病态的问题。病态问题通常需要预处理、或通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）、或者重新寻找出现病态的原因，改变原问题的提法。

问题提出：考虑一个高次的代数多项式

∏=-=---=201)()20)...(2)(1()(k k x x x x x p

（E1-1）

显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简单的）。考虑该多项式方程的一个扰动

0)(19=+x x p ε

(E1-2) 其中ε是一个非常小的数。这相当于是对（E1-1）中19

x 的系数作一个小的扰动。我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。 实验内容：为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数

u ＝roots （a ）

其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程 0...1121=++++-n n n n a x a x a x a

的全部根，而函数b=poly(v)的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。

ve=zeros(1,21);

ve(2)=0.01;（将ve 中的第二个分量改为0.01）

roots(poly(1:20))+ve)

上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“0.01”即是（E1-2）中的ε。

实验要求：

（1）选择充分小的ε，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数ε很小，我们自然感觉（E1-1）和(E1-2)的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？

（2）将方程（E1-2）中的扰动项改成18

x ε或其他形式，实验中又有怎样的现象出现？ 实验1.2 误差传播与算法稳定性

问题提出：考虑一个简单的由积分定义的序列 ...2,1,101==?-n dx e x I x n n

（E1-3） 显然,...2,1,0=>n I n .当1=n 时，e dx xe I x /11011==

?-。而对于2≥n 时，利用分部积分易得

,...3,2,1|11

01110110

1=-=-==-----??n nI dx e nx e x dx e x I n x n x n x n n 另一方面，我们有

1111100011(1)1n n x n n x e dx I x e dx x dx e n n --≤≤=≤=++???. （E1-4）

实验内容：由以上递推关系，我们可得到计算序列{n I }的两种方法。

（I ）：,...3,2,1,/111=-==-n nI I e I n n

（E1-5） （II ）：2,3,...2,1,,1,01--=-==-N N N n n

E E E n n N （E1-6） 实验要求：

（1）分别用算法（I ）、（II ）并在计算中分别采用5位、6位和7位有效数字，请判断哪种算法能给出更精确的结果。

（2）两种算法的优劣，与你的第一感觉是否吻合？解释实验的结果。算法（I ）中的1I 的计算误差为1e ，由1I 递推计算n I 的误差为n e ；算法（II ）中N I 的计算误差为N ε，由N I 递推计算)(N n I n <的误差为n ε。如果在上述两算法中都假定后面的计算不再引入其他误差， 试给出n e 与1e 的关系和n ε与N ε的关系。

数值实验综述：非线性方程的数值解法是科学与工程计算的一个重要课题，在实际问题中有广泛的应用，特别在各种非线性问题的科学计算中更显出它的重要性。本章讨论非线性方程的数值解法，主要讲授了解非线性方程的二分法、迭代法、迭代加速、Newton 法割线法和抛物线法等。本章数值实验的中心就是通过具体的数值实验来研究算法的可行性，如算法的收敛性，初值的选择，收敛速度以及算法的复杂性与稳定性等。

实验2.1 迭代法、初始值与收敛性

初步认识非线性问题的迭代法与线性问题迭代法的差别，探讨迭代法及初始值与迭代收敛性的关系。

问题提出：迭代法是求解非线性方程及方程组的基本思想方法，与线性方程的情况一样，其构造方法可以有多种多样，但关键是怎样才能使迭代收敛且有较快的收敛速度。

实验内容：考虑一个简单的代数方程

012=--x x ，

针对上述方程，可以构造多种迭代法，如121-=+n n x x ,n

n x x 111+=+,11+=+n n x x 等。在实轴上取初始值0x ，分别用以上迭代作实验，记录各算法的迭代过程。

实验要求：

（1）取定某个初始值，按如上迭代格式进行计算，它们的收敛性如何？重复选取不同的初始值，反复实验。请读者自行设计一种比较形象的记录方式（如利用Matlab 的图形功能），分析三种迭代法的收敛性与初值选取的关系。

（2）对三个迭代法中的某一个，取不同的初始值进行迭代，结果如何？试分析迭代法对不同的初值是否有差异？

（3）线性方程组迭代法的收敛性是不依赖初始值选取的。比较线性与非线性问题迭代的差异，有何结论和问题。

实验2.2 算法设计与比较

非线性方程的数值求解方法很多，基本的思想是线性化。不同的方法效果如何，要靠计算的实践来分析、比较。

实验内容：考虑算法：（l ）牛顿法；（2）割线法

分别编写它们的Matlab 程序。

实验要求：（1）用上述各种方法，分别计算下面的两个例子。在达到精度相同的前提下，比较其迭代次数、浮点运算次数和CPU 时间等。

（I ）3231080,00x x x +-==取

（II ）22123081(0.1)sin 1.060,

0x x x x -+++==取。 （2）取其它的初值0x ，结果如何？反复选取不同的初值，比较其结果。

（3）总结归纳你的实验结果，试说明各种方法的特点。

数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。求解方法大致可分为直接法和迭代法两大类。直接法指在不考虑舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss 消去法是目前最基本和常用的方法。迭代法是逐次逼近的方法，可用于求解大型特别是大型稀疏的方程组。迭代法的优点是可以利用系数矩阵的特殊结构，尽可能地减少计算量与存储量，但需要判断收敛性。

实验3.1 Gauss 消去法与LU 分解

Gauss 消去法等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。对一般的正定对称矩阵，采用平方根方法可不选主元。

实验内容：请用不选主元的LU 分解和列主元LU 分解求解下列线性方程组 Ax=b 。

实验要求：

（1）比较这两种方法的结果；

（2）增加方程阶数，再比较分析所得结果。

实验3.2 直接法的时间复杂性试验

不同的算法所需要的计算量有所不同，因此所花费的时间也不同。解n 阶线性方程组的直接法Gauss 消去法或者LU 分解法所花计算量大约为：32

3n ，而利用三角分解求逆所花工作量大约为：383n 。故利用逆求解线性方程组远比Gauss 消去法或LU 分解更费时。如果是

三对角方程组，则用追赶法所花工作量约为：87n -。应当注意的是，随着计算机的运算速度越来越快，纯粹计算所花时间比重下降，存取和通讯的量变得越来越重要，故存储量和存取方式也是影响计算速度的一个重要因素。

实验内容：请分别用不同方法求解实验3.1的线性方程组 Ax=b ，记录各方法所用时间，并与各方法的不同工作量进行比较。

实验要求：

（1）分别用\x A b =，()*x inv A b =计算，并分别记录所花费的CPU 时间，进行分析比较。

（2）取300=n ，随机生成A 的一条主对角线元素和二条次对角线元素，使A 为严格对角占优的三对角阵；随机生成b ；其中三条对角线元素分别用三个一维数组存储；用Matlab 语言自编M 文件分别用\x A b =，()*x inv A b =和追赶法计算这个三对角方程组；并分别记录所花费的CPU 时间；分析结果，得出你的结论。

实验3.3 线性方程组的迭代法

线性方程组的迭代法是一类重要方法，如Jacobi 迭代法、Gauss-Seidel 迭代法、SOR 方法等都是一些经典方法，虽已不是最优的方法，但还有其基础作用。本实验通过具体的数值例子来研究算法的收敛性，初值的选择，收敛速度等。

实验内容3.3：线性方程组的迭代法比较

实验要求：

（1）编写Jacobi 迭代法、Gauss-Seidel 迭代法、SOR 迭代法的M 文件。

（2）给定2020?∈R A 为五对角矩阵

???????????????

?????????????????------------------321

412132141412132141412132141412132141213 选取不同的初始向量)0(x 及右端面向量b ，给定迭代误差要求，利用已编写的Jacobi 迭代法和GS 迭代法程序求解，观察得到的序列是否收敛？若收敛，通过迭代次数分析计算结果并得出你的结论。用已编写的SOR 迭代法程序，对于所选取的初始向量)0(x 及右端项b 进行求

解，松驰系数ω取1<ω<2的不同值，在5)1()(10-+≤-k k x

x 时停止迭代，通过迭代次数分析计算结果并得出你的结论。

数值实四

插值与拟合是数值逼近的两种主要形式，数值实验的目的包括以下几个方面：（1）各种插值多项式的构造，高阶插值的Runge 现象，分段低次插值特别是样条插值的实现。 (2) 插值与拟合的区别与联系。

t x

-2 实验内容：(1) 分别用Lagrange 插值法和Newton 插值法编程求解。

(2) 求出插值多项式的系数，对比计算结果。

实验要求：（1）分析对比这两种插值方法的计算过程与结果。

（

2）绘制插值函数曲线，并分析数据点的变化趋势。

实验4.2 三次样条插值

（2）计算出各插值节点的弯矩值。

实验要求：绘制函数()f x 曲线，并与分段线性插值效果进行比较。

实验4.3龙格（Runge ）现象

考虑在一个固定的区间上高次等距插值多项式逼近一个函数

22511)(x

x f += 实验内容：考虑区间]1,1[-的一个等距划分，分点为

n i n

i x i ,...,2,1,0,21=+

-= 则拉格朗日插值多项式为 ∑=+=n

i i

i

n x l x x L 02)(2511)( 其中的)(x l i ，n i ,...,2,1,0=是n 次拉格朗日插值基函数。

实验要求：（l ）选择不断增大的分点数目,...3,2=n ，画出原函数)(x f 及插值多项式函数)(x L n 在]1,1[-上的图像，比较并分析实验结果。

（2）选择其他的函数，例如定义在区间[-5，5]上的函数

x x g x

x x h arctan )(,1)(4=+= 重复上述的实验看其结果如何。

（3）区间],[b a 上切比雪夫点的定义为

1,...,2,1,)1(2()12(cos 22+=???

? ??+--++=n k n k a b a b x k π 以121,...,,+n x x x 为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果。

实验4.4 曲线逼近方法的比较

曲线的拟合和插值，是函数逼近的基本方法，每种方法具有各自的特点和特定的适用范围，实际工作中合理选择方法是重要的。

实验要求：（1）分别画出三种方法的近似曲线；（2）分析比较拟合的结果与拉格朗日插值及样条插值的结果。

数值微积分是科学与工程计算中常见内容，内容丰富，本章只介绍了其中最基本的部分。如高精度的Gauss 积分，奇异积分，多元函数积分，积分方程数值解等常见内容并没有介绍。

实验5.1 多种数值求积公式计算定积分的比较

实验内容：计算积分4

2013()exp( 1.5)x x x dx --?的近似值

实验要求：（1）用复化梯形公式，（2）自动选取步长的梯形法，（3）Romberg 算法计算，要求绝对误差限为810ε-=

实验5.2 奇异积分探索

实验内容：考虑积分1

01sin x dx x ?. 实验要求：（1）用MATLAB 函数explot 在11x -≤≤上画被积函数的图像；

（2）用MATLAB 的符号计算工具箱，精确计算该积分；

（3）如果直接用MATLAB 函数quad 和quadl 数值计算该积分有什么问题？如何设法克服?

微分方程是科学与工程计算中最常见的形式之一。微分方程的解也非常复杂，尽管有很多数值方法，但数值方法也常常遇到很多困难，如稳定性、收敛性、计算量和计算速度、刚性等问题。常微分方程初值问题的数值解法又是边值问题偏微分方程数值解的基础。本章重点训练初值问题数值解的各种方法和稳定性等问题。

实验6.1 数值方法的绝对稳定性

通过不同步长的选择体会各种不同方法稳定域的概念。

实验内容：

对于微分方程的初值问题：()'50y x y =-，[0,]x b ∈，()01y =

分别采用向前Euler 法、向后Euler 法、改进Euler 法和经典四级四阶Runge-Kutta 法编程求数值解。

实验要求：（1）步长分别取0.1,0.05,0.025,0.001h =；

（2）区间长度分别取1,10,100b =

（3）与MATLAB 函数ode45和ode23的计算结果以及精确解对比。

实验6.2 求解刚性问题

刚性问题是常微分方程组中的一个“病态问题”，即它反映的是稳定性要求必须小步长，但求解范围是大区间之间的一对矛盾。

实验内容：用改进Euler 方法求解下列微分方程组的初值问题：

11212122

'1000.25999.750.5,1,'999.751000.250.5,1y y y y y y y y =-++=??=-+=-?，050x << 实验要求：

（1）编写改进Eule 解微分方程组的MATLAB 程序，实验得到最小稳定步长；

（2）编写向后Euler 法解此问题的MATLAB 程序，实验得到最小稳定步长；

（3）比较二者效率，并与MATLAB 函数oda15s 求解效率作比较。