数值分析上机作业1-1教案.doc

1、 1 数值计算方法上机题目 1 1、 实验 1. 病态问题 实验目的 ： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出 ： 考虑一个高次的代数多项式 201 )()20) . . . (2)(1()( k kxxxxxp（ E1-1） 显然该多项式的全部根为 l

2、， 2， 20，共计 20 个，且每个根都是单重的（也称为简单的）。现考虑该多项式方程的一个扰动 0)( 19 xxp (E1-2) 其中 是一个非常小的数。这相当于是对（ E1-1）中 19x 的系 数作一个小的扰动。我们希望比较（ E1-1）和（ E1-2）根的差别，从而分析方程（ E1-1）的解对扰动的敏感性。 实验内容 ： 为了实现方便，我们先介绍两个 Matlab 函数：“ roots”和“ poly”，输入函数 u roots （ a ） 其中若变量 a 存储 1n 维的向量，则该函数的输出 u 为一个 n 维的向量。设 a 的元素依次为 121 ,., naaa ，则输出 u 的

3、各分量是多项式方程 0. . . 1121 nnnn axaxaxa 的全部根，而函数 b=poly(v) 的输出 b 是一个 n 1 维变量，它是以 n 维变量 v 的各分量为根的多项式的系数。可见“ roots”和“ Poly”是两个互逆的运算函数 . ve=zeros(1,21); ve(2)=ess; roots(poly(1:20)+ve) 上述简单的 Matlab 程序便得到（ E1-2）的全部根，程序中的“ ess”即是（ E1-2）中的 。 实验要求 ： （ 1）选择充分小的 ess，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数 很小，我们自然感觉（ E1-1）

4、和 (E1-2)的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？ （ 2）将方程（ E1-2）中的扰动项改成 18x 或其他形式，实验中又有怎样的现象出现？ 2 实验步骤： （ 1）程序 function t_charpt1_1 clc result=inputdlg(请输入扰动项 :在 0 20之间的整数 :,charpt 1_1,1,19); Numb=str2num(char(result); if(Numb20)|(Numb0)errordlg(请输入正确的扰动项 :0 20之间的整数 !);return;end result=inputdlg(请输

5、入 (0 1)之间的扰动常数 :,charpt 1_1,1,0.00001); ess=str2num(char(result); ve=zeros(1,21); ve(21-Numb)=ess; root=roots(poly(1:20)+ve); x0=real(root); y0=imag(root); plot(x0,y0, *); disp(对扰动项 ,num2str(Numb),加扰动 ,num2str(ess),得到的全部根为 :); disp(num2str(root); 二、 实验结果分析 ess 分别为 1e-6,1e-8.1e-10,1e-12. 对扰动项 19 加扰动

6、1e-006 得到的全 部根为 : 21.3025+1.56717i 21.3025-1.56717i 18.5028+3.6004i 18.5028-3.6004i 15.1651+3.76125i 15.1651-3.76125i 12.4866+2.88278i 12.4866-2.88278i 10.5225+1.71959i 10.5225-1.71959i 9.04485+0.594589i 9.04485-0.594589i 7.9489+0i 7.00247+0i 5.99995+0i 5+0i 4+0i 3+0i 2+0i 1+0i 对扰动项 19 加扰动 1e-010 得到

7、的全部根为 : 19.9953+0i 19.0323+0i 17.8696+0i 17.2186+0i 15.4988+0.0211828i 15.4988-0.0211828i 13.7707+0i 13.1598+0i 11.9343+0i 11.029+0i 9.99073+0i 9.00247+0i 7.99952+0i 7.00007+0i 5.99999+0i 5+0i 4+0i 3+0i 2+0i 1+0i ess 分别为 1e-6,1e-8.1e-10,1e-12 的图像如下： 3 从实验的图形中可以看出，当 ess 充分小时，方程 E.1.1 和方程 E.1.2 的解相差很小

8、，当ess 逐渐增大时，方程的解就出现了病态解，这些解都呈现复共轭性质。 （ 2） 将扰动项加到 x18上后， ess=1e-009 时方程的解都比较准确，没有出现复共轭现象。ess=1e-008 时误差与 x19(ess=1e-009)时相当，即扰动加到 x18上比加到 x19小一个数量级。对x8的扰动 ess=1000 时没有出现复共轭，误差很小；对 x 的扰动 ess=10e10 时没有出现复共轭，误差很小。因此，扰动作用到 xn上时， n 越小，扰动引起的误差越小。 2、 实验 2。 多项式插值的振荡现象，即插值的龙格（ Runge）现象 问题提出 ： 考虑在一个固定的区间上用插值逼近

9、一个函数。显然，拉格朗日插值中使用的节点越多，插值多项式的次数就越高、自然关心插值多项式的次数增加时， )(xLn 是否也更加靠近被逼近的函数。龙格给出的一个例子是极著名并富有启发性的。设区间 1,1 上函数 2251 1)( xxf 实验内容 ： 考虑区间 1,1 的一个等距划分，分点为 ninix i ,.,2,1,0,21 则拉格朗日插值多项式为 ni iin xlxxL 0 2 )(251 1)( 其中的 )(xli ， ni ,.,2,1,0 是 n 次拉格朗日插值基函数。 实验要求 ： （ l）选择不断增大的分点数目 ,.3,2n ，画出原函数 )(xf 及插值多项式函数 )(xL

10、n 在1,1 上的图像，比较并分析实验结果。 （ 2）选择其他的函数，例如定义在区间 -5， 5上的函数 4 xxgxxxh a r c ta n)(,1)( 4 重复上述的实验看其结果如何。 （ 3）区间 , ba 上切 比雪夫点的定义为 1, . . . ,2,1,)1(2( )12(c o s22 nknkababx k 以 121 ,., nxxx 为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果。 实验步骤： （ 1） 试验程序： function y=Lagrange(x0, y0, x); % Lagrange插值 n= length(x0); m=length(x);

11、for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if(j = k) p = p*(z - x0(j)/(x0(k) - x0(j); end end s = s + p*y0(k); end y(i) = s; end function t_charpt2 promps = 请选择实验函数，若选 f(x),请输入 f,若选 h(x),请输入 h,若选 g(x),请输入 g:; titles = charpt_2; result = inputdlg(promps,charpt 2,1,f); Nb_f = char(result); if

12、(Nb_f = f return;end result = inputdlg(请输入插值结点数 N:,charpt_2,1,10); Nd = str2num(char(result); if(Nd 1)errordlg(结点输入错误！ );return;end switch Nb_f 5 case f f=inline(1./(1+25*x.2); a = -1;b = 1; case h f=inline(x./(1+x.4); a = -5; b = 5; case g f=inline(atan(x); a = -5; b= 5; end x0 = linspace(a, b, Nd+

13、1); y0 = feval(f, x0); x = a:0.1:b; y = Lagrange(x0, y0, x); fplot(f, a b, co); hold on; plot(x, y, b-); xlabel(x); ylabel(y = f(x) o and y = Ln(x)-); 增大分点 n=2， 3，时，拉格朗日插值函数曲线如图所示。 n=3 n=6 n=7 n=8 从图中可以看出，随着 n 的增大，拉格朗日插值函数在 x=0 附近较好地逼近了原来的函数 f(x)，但是却在两端 x= -1 和 x=1 处出现了很大的振荡现象。 通过 分析图形，可以看出，当 n 为奇数时

14、，虽然有振荡，但振荡的幅度不算太大， n 为偶数时，其振荡幅度变得很大。 （ 2） 将原来的 f(x)换为其他函数如 h(x)、 g(x)，结果如图所示。其中 h(x), g(x)均定义在 -5，5区间 上， h(x)=x/(1+x4)， g(x)=arctan x。 6 h(x), n=7 h(x), n=8 h(x), n=9 h(x), n=10 g(x), n=8 g(x), n=9 g(x), n=12 g(x), n=13 分析两个函数的插值图形，可以看出： 随着 n 的增大，拉格朗 日插值函数在 x=0 附近较好地逼近了原来的函数 f(x)，但是却在两端 x= -5 和 x=5

15、处出现了很大的振荡现象。 通过图形 可以看出，当 n 为偶数时，虽然有振荡，但振荡的幅度不算太大， n 为奇数时，其振荡幅度变得很大。原因和上面 f(x)的插值类似， h(x)、 g(x)本身是奇函数，如果 n 为偶数，那么 Lagrange 插值函数 Ln(x)的最高次项 xn-1是奇次幂，比较符合 h(x)、 g(x)本身是奇函数的性质；如果 n 为奇数，那么 Lagrange 插值函数 Ln(x)的最高次项 xn-1是偶次幂，与 h(x)、 g(x)本身是奇函数的性质 相反，因此振荡可能更7 剧烈。 3、 实验 3。 样条插值的收敛性 问题提出 ： 一般的多项式插值不能保证收敛性，即插值

16、的节点多，效果不一定就好。对样条函数插值又如何呢？理论上证明样条插值的收敛性是比较困难的，也超出了本课程的内容。通过本实验可以验证这一理论结果。 实验内容 ： 请按一定的规则分别选择等距或者非等距的插值节点，并不断增加插值节点的个数。考虑实验 2.中的函数或选择其它你有兴趣的函数，可以用 Matab 的函数 “spline”作此函数的三次样条插值。在较新版本的 Matlab 中，还提供有 spline 工具箱，你 可以找到极丰富的样条工具，包括 B-样条。 实验要求 ： （ 1）随节点个数增加，比较被逼近函数和样条插值函数误差变化情况。分析所得结果并与拉格朗目多项式插值比较。 （ 2）样条插值

17、的思想最早产生于工业部门。作为工业应用的例子，考虑如下问题：某汽车制造商用三次样条插值设计车门的曲线，其中一段的数据如下： kx 0 1 2 3 4 5 6 7 8 9 10 ky 0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 ky 0.8 0.2 实验步骤： （ 1）程序： function t_charpt2 promps = 请选择实验函数，若选 f(x),请输入 f,若选 h(x),请输入 h,若选 g(x),请输入 g:; titles = charpt_2; result = inputdlg(promps,charpt

18、 2,1,f); Nb_f = char(result); if(Nb_f = f return;end result = inputdlg(请输入插值结点数 N:,charpt_2,1,10); Nd = str2num(char(result); if(Nd 1)errordlg(结点输入错误！ );return;end switch Nb_f case f f=inline(1./(1+25*x.2); a = -1;b = 1; case h f=inline(x./(1+x.4); a = -5; b = 5; case g 8 f=inline(atan(x); a = -5; b

19、= 5; end x0 = linspace(a, b, Nd+1); y0 = feval(f, x0); x = a:0.1:b; cs = spline(x0, y0); y = ppval(cs, x); plot(x0, y0, o); hold on; plot(x, y, k-); xlabel(x); ylabel(y = f(x) o and y = Spline(x)-); 实验结果： 如图所示。 f(x), n=5 n=10 n=20 h(x), n=5 h(x), n=10 n=20 g(x), n=5 n=10 n=20 图中可以看出，由于其采用了分段三次多项式拟合的方法，随着三次样条插值 的 插值结点的增加，并没有出现振荡现象。 （ 2） 程序： x0=0:10; y0=0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29; x=0:0.1:10; pp=csape(x0,y0,complete,0.8 0.2); y = ppval(pp, x); plot(x0, y0, o); hold on; plot(x, y, k-); xlabel(x); ylabel(y = f(x) o and y = Spline(x)-); 9 车门的曲线 如下图：