数值分析课程设计报告.pdf

《数值分析课程设计报告.pdf》由会员分享，可在线阅读，更多相关《数值分析课程设计报告.pdf（20页珍藏版）》请在三一文库上搜索。

1、1 数 值 分 析 课 程 设 计 报 告 摘要（中）：本文建立在数值分析的理论基础上，能够在Matlab 环境中运行， 给出了理论分析、程序清单以及计算结果。更重要的是，还有详细的对算法的 框图说明。首先运用Romberg积分方法对给出定积分进行积分，然后对得到的 结果用插值方法，分别求出Lagrange 插值多项式和 Newton 插值多项式，再运 用最小二乘法的思想求出拟合多项式，最后对这些不同类型多项式进行比较， 找出它们各自的优劣。 Summary(Enlish):This essay is based on numerical analysis ,and could be oper

2、ated in Matlab environment, including theroy analysis, program and results. What s more, there are detailed diagram which shows how the algorithm works. It first uses the integrating method of Romberg , which is an improved trapezoidal integration, to solve the given definite integral,then we create

3、 Lagrange s interpolation polynomial and Newton s interpolation polynomial. And according to least square method, curve fitting polynomial is created. At the last part of the essay, we compare these different patterns of polynomial, founding their distinctive advantages and disadvantages. 主题词：Romber

4、g积分，插值方法， Langrange 插值多项式， Newton插值多项式， 拟合多项式。 2 一. 问题提出： 已知椭圆的周长可以表示成s=a 2 0 22 cos1d(00.0001 D=zeros(n,n); D(:,1)=y; for j=2:n for k=j:n D(k,j)=(D(k,j-1)-D(k-1,j-1)/(x(k)-x(k-j+1); end end D C=D(n,n); for k=(n-1):-1:1 C=conv(C,poly(x(k); m=length(C); C(m)=C(m)+D(k,k); end 构造 n阶零矩阵 D 把 y 的值存入 D 的第一

5、列 依次计算第k 阶均差，保留在D 的第 (k+1)列 根据 (*) 式利用差商即D 对角线上值元素求得Newton 插值多项式 8 运行结果： D = Columns 1 through 7 6.2989 0 0 0 0 0 0 6.3456 0.4670 0 0 0 0 0 6.4223 0.7670 1.5000 0 0 0 0 6.5274 1.0510 1.4200 -0.2667 0 0 0 6.6592 1.3180 1.3350 -0.2833 -0.0417 0 0 6.8153 1.5610 1.2150 -0.4000 -0.2917 -0.5000 0 6.9936 1

6、.7830 1.1100 -0.3500 0.1250 0.8333 2.2222 7.1919 1.9830 1.0000 -0.3667 -0.0417 -0.3333 -1.9444 7.4082 2.1630 0.9000 -0.3333 0.0833 0.2500 0.9722 Columns 8 through 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -5.9524 0 4.1667 12.6488 ans = Columns 1 through 7 12.6488 -51.4881 87.9514 -81.7083 44.6592 -14.8231 4.394

7、6 Columns 8 through 9 -0.2795 6.2940 将结果整理，得到Newton 插值多项式为： y=12.6488x 8-51.4881x7+87.9514x6 -81.7083x 5+44.6592x4-14.8231x3+4.3946x2-0.2795x+ 6.2940 运用 Matlab 画出 Newton 插值多项式的拟合效果图，用叉“x”表示原数据点， 编写 Matlab 程序 newpolyplot.m function newpolyplot x=0.1:0.1:0.9; y=6.2989 6.3456 6.4223 6.5274 6.6592 6.815

8、3 6.9936 7.1919 7.4082; c,d=newpoly(x,y); y0=polyval(c,x); plot(x,y,x,x,y0) 9 可以看出，两个插值方法得出的结果相同，不过，Newton 插值要比Lagrange 插值好， 由于 Lagrange 插值的n 次插值基函数lk(x)(k=0,1,2, n)都依赖于全部插值结点，利用 公式很容易得到插值多项式，但在增加或减少结点时，插值基函数lk(x)(k=0,1,2, n)也 随之变化，必须全部重新计算。而克服这个缺点的有效方法之一就是构造Newton 插值 那样的均差表，在插值点变化时，有些数据可以保留。 ?拟合多项式

9、 曲线拟合问题与函数插值问题不同，它不要求曲线通过所有已知点，只要求得到的 近似函数能反映数据的基本关系。因此，曲线拟合的过程比插值过程得到的结果更能反 映客观实际。 在某种意义上， 曲线拟合更具有实用价值，因为实际问题中所提供的数据 往往很多，如果用插值法势必要得到次数很高的插值多项式，导致计算上的很多麻烦。 一般总是希望各观察数据与拟合曲线的偏差的平方和最小，这样就能使拟合曲线更 接近于真实函数。这个原理就称为最小二乘原理。用最小二乘原理作为衡量“曲线拟合 优劣”的准则称为曲线的最小二乘法。 最小二乘法的思想是： 设有 N 个数据点 （xk,yk） ， 并给定 M 个线性独立函数f j(x

10、) 。 为求 M 个稀疏，使用由线性组合形成的函数f(x) ，表示为： M j jj xfcxf 1 求解最小误差平方和，表示为： N k k M j kjj M k kkM yxfcyxfcccE 1 2 11 2 21 ,., 为求解E 的最小值，每个偏导数必须为零（即MicE i ,.,2, 1,0/） ，这样 可得到如下方程组： 10 Mixfyxfc N k kik M j kjj ,.,2, 1,0 11 交换上式中的求和顺序，可得一个M*M线性方程组，未知数是系数cj，称为正 规方程组： M j N k kkij N k kjkj Miyxfcxfxf 111 ,.2, 1,*

11、将上面的正规方程组表示成矩阵形式可减少不必要的计算量。 方程组 * 的系数矩阵是一个对称矩阵，并且是正定的，可以唯一地确定系数cj， 根据这些特点，构造如下矩阵： f1(x1) f2(x1) ,fM(x1) f1(x2) f2(x2) ,fM(x2) F= f1(x3) f2(x3) ,fM(x3) , f1(xN) f2(xN) ,fM(xN) f1(x1) f1(x2) ,f1(xN) f2(x1) f2(x2) ,f2(xN) F = f3(x1) f3(x2) ,f3(xN) , fM(x1) fM(x2) ,fMxN) 乘积 F F 是一个 M*M矩阵，方程组 * 转化为求解线性方程

12、组： YFFCF （求解系数矩阵C） 程序框图： 编写 Matlab 程序 lspoly.m: function c=lspoly(x,y,M) n=length(x); B=zeros(1:M+1); F=zeros(n,M+1); 构造 M+1 维零向量B，N 行 M+1 列的零矩阵F 将 x k-1 依次存放到F 的第 k 列中 求解正规方程组F FC=F Y 11 for k=1:M+1 F(:,k)=x.(k-1); end A=F*F; B=F*y; c=(AB); 运行结果为： x=0.1:0.1:0.9; y=6.2989 6.3456 6.4223 6.5274 6.6592

13、 6.8153 6.9936 7.1919 7.4082; lspoly(x,y,2) ans = 6.2624 0.1864 1.2147 lspoly(x,y,4) ans = 6.2844 -0.0200 1.6838 -0.2704 -0.0382 lspoly(x,y,6) ans = 6.2829 0.0069 1.5243 0.1506 -0.5697 0.3013 -0.0556 lspoly(x,y,8) ans = Columns 1 through 8 6.2940 -0.2795 4.3940 -14.8201 44.6498 -81.6911 87.9326 -51.

14、4769 Column 9 12.6460 整理得到偶次拟合多项式： 2 阶拟合多项式：y=6.2624+0.1864x+1.2147x 2 4 阶拟合多项式：y=6.2844-0.02x+1.6838x 2-0.2704x3-0.0382x4 6 阶拟合多项式：y=6.2829+0.0069x+1.5243x 2+0.1506x3-0.5697x4+0.3013x5-0.0556x6 8 阶拟合多项式： y=6.2940-0.2795x+4.3940x 2-14.8201x3 +44.6498x 4-81.6911x5+87.9326x6-51.4769x7 +12.6460x 8 为检验每

15、个拟合多项式的效果，编写程序lspoly2.m： function Y=lspoly2(M) x=0.1:0.1:0.9; y=6.2989 6.3456 6.4223 6.5274 6.6592 6.8153 6.9936 7.1919 7.4082; C=(lspoly(x,y,M); %存放拟合多项式的系数 n=length(x); X=zeros(n,M+1); for k=1:M+1 X(:,k)=x.(k-1); end Y=(X*C); 12 将结果列出如表： 实际值6.2989 6.3456 6.4223 6.5274 6.6592 6.8153 6.9936 7.1919 7

16、.4082 2 阶6.2932 6.3483 6.4277 6.5313 6.6593 6.8115 6.9881 7.1889 7.4141 4 阶6.2989 6.3455 6.4223 6.5275 6.6592 6.8152 6.9936 7.1920 7.4082 6 阶6.2989 6.3456 6.4223 6.5274 6.6592 6.8153 6.9936 7.1919 7.4082 8 阶6.2989 6.3456 6.4223 6.5274 6.6592 6.8153 6.9936 7.1919 7.4082 显然， 6 阶拟合多项式的拟合结果已与实际值吻合了。 既然到

17、 6 阶拟合时，拟合结果已经令人满意就取6 阶拟合多项式： y=6.2829+0.0069x+1.5243x 2+0.1506x3-0.5697x4+0.3013x5-0.0556x6 三结果与分析 : 第?小题的运算过程清楚地表明：Romberg 积分方法高速有效，易于编制程序， 适合于计算机上操作。根据 Romberg积分方法构造出来的序列nmTm , 1 , 0,，其收 敛速度比变步长求积法更快。这是因为，在Romberg 积分方法中，采用了提高阶数与 减少步长这两种提高精度的措施。但它有一个明显的缺点是：每当把区间对分后，就 要对被积函数f(x)计算它在新分点处的值，而这些函数值的个数

18、是成倍增加的。 第?小题先后使用了Lagrange 插值方法和Newton 插值方法， Lagrange 插值多项 式具有明显地对称性，易于记忆。 然而， 正是由于这种对称性使得公式中的每一项与所 有的插值结点有关。因此，如果增加一个插值结点，则Lagrange 插值多项式中的每一 项多需要改变。由Newton 插值公式可以看出，当增加一个插值结点时，前面已有的各 项均不变，只需在后面增加一项即可。并且，Newton 插值多项式中，每一项的系数就 是各阶差商，它减少了不必要的重复计算，比Lagrange 插值多项式计算量小，且便于 程序设计。 在第 ?小题中，得到Lagrange 插值多项式和

19、Newton 插值多项式的结果是一样的。 这是因为通过相同结点的插值多项式是唯一的，因此，Lagrange 插值多项式与Newton 插值多项式本质上是同一个多项式，只不过形式不同而已，这恰验证了代数插值问题的 唯一有解性： 不管构造插值多项式采用哪种方式，只要满足给定的n+1 个结点上的函数 yi=f(xi),i=0,1,n.将插值多项式化成一般多项式Pn(x)=anx n+a n-1x n-1+, +a 0后，其系数是 唯一的。 第 ?小题所讨论的插值问题，是用一个多项式来近似代替列表函数，并且要求多项 式通过列表汉书中给定的数据点。而在实际应用中，还经常会遇到这样的问题，为解析 式子比较

20、复杂的函数寻找一个多项式来近似代替它，并要求其误差在某种度量意义下最 小。 另外，插值法虽然在一定程度上可以解决根据给定函数值求函数的近似表达式问题， 但同时亦存在着明显的缺陷。首先， 因为由实验提供的数据往往带有测试误差，个别数 据的误差还可能很大，如果要求近似曲线严格地经过所有给定的数据点，就会使曲线保 留这些误差， 从而失去原数据表示的规律；其次，即使插值多项式通过了给定的数据点， 在这些给定数据点上的误差很小，但不能保证在其他点上的误差也很小，所以， 在两个 插值结点之间，插值多项式Pn(x)并不一定能很好地逼近被插值函数 f(x)，有时差异甚 至很大。同时，当实验数据很多时，用插值法

21、得到近似表达式明显地缺乏实用价值。 第 ?小题中， 取 6 阶拟合多项式比较好，因为高阶多项式的拟合会引起数值的不稳 定，使用最小二乘法多项式拟合非线性数据的方法固然具有吸引力，但如果数据不具有 多项式特性， 则求出的曲线可能产生大的振荡。这种现象在高阶多项式情况下更容易发 生，由于这个原因，一般很少使用超过6 阶的多项式。正因为这样，所以才会选6阶拟 合多项式。 13 参考文献： John H.Mathews,Kurtis D.Fink 著，陈渝等译。Numerical Methods Using MATLAB (Third Edition) 北京：电子工业出版社. 张池平主编 .计算方法

22、.哈尔滨：哈尔滨工业大学出版社. 王兵团，桂文豪著.数学实验基础 .北京：北方交通大学出版社. 杨万利等著 .数值分析教程.北京：国防工业出版社. 康晓东主编 .数值算法与非数值算法. 北京：电子工业出版社 徐士良编著 .数值分析与算法.北京：机械工业出版社. 林成森编著 .数值计算方法.北京：科学出版社. 张池平，施云慧著.计算方法 . 北京：科学出版社. 西安郵電學院 数值分析课程设计报告书 系 部 名 称：应用数理系 学 生 姓 名： 14 专 业 名 称：信息与计算科学 班级：信息 070x 时间：2 15 实验 X 19 次方系数对方程根的影响 一、问题的提出 考虑代数多项式 20

23、1 ( )(1)(2)(20)() k P xxxxxk, 显然它的全部根 为 1,2, , ,20; 如果 19 次方的系数有扰动 , 那么对方程的解有何影响. 二、实验内容 直接利用 Matlab 中的 roots()和 poly() 函数, 选择不同的扰动做数值计算. 主要步骤 :ve=zeros(1,21);ve(2)=ess;roots(poly(1:20)+ve). Matlab 程序: 见附件 bingtaiwenti.m 三、实验结果 1. 取扰动分别为 369 10 ,10 ,10所得结果如下 : 3 10 6 10 9 10 27.0817+5.03812i 21.3025

24、+1.56717i 19.952+0i 27.0817-5.03812i 21.3025-1.56717i 19.2293+0i 19.5337+9.1664i 18.5028+3.6004i 17.6573+0.692896i 19.5337-9.1664i 18.5028-3.6004i 17.6573-0.692896i 13.8235+7.77167i 15.1651+3.76125i 15.4524+0.875524i 13.8235-7.77167i 15.1651-3.76125i 15.4524-0.875524i 10.7211+5.4609i 12.4866+2.88278

25、i 13.3527+0.486992i 10.7211-5.4609i 12.4866-2.88278i 13.3527-0.486992i 8.91282+3.47317i 10.5225+1.71958i 11.8578+0i 8.91282-3.47317i 10.5225-1.71958i 11.0427+0i 7.69268+1.89884i 9.04487+0.59455i 9.9916+0i 7.69268-1.89884i 9.04487-0.59455i 9.00201+0i 6.75761+0.654714i 7.94891+0i 7.99952+0i 6.75761-0.

26、654714i 7.00247+0i 7.00009+0i 5.95208+0i 5.99995+0i 5.99999+0i 5.00061+0i 5+0i 5+0i 4+0i 4+0i 4+0i 3+0i 3+0i 3+0i 2+0i 2+0i 2+0i 1+0i 1+0i 1+0i 16 2. 绘制相应的图形 , 即将以上结果可视化 . 图 2. 6 10 图 1. 3 10 17 四、实验结果分析 观测现象 : 由结果可以观察到误差是不可避免的, 并且扰动的减小不能使得所有根除的 误差减小 , 还可以观察到16x处的误差总是最大 , 而靠近1x处的误差较小 . 误差分析 : 原问题数学上

27、描述为 : 19 ( )0P xx的求解问题 , 也可等价为 : 2019 ( ,)0P xxx 那么要考虑的问题是对方程的某个解x 的影响 , 我们不妨将 x 看成的函数 , 从而可以求得此问题的条件数为: 18 2020 11 ( ,) () p ij ji P dxax Cx P x dx xj x 可以求得 20 个根处的条件数为 : 图 3. 9 10 18 19次方系数变化对根的影响问题条件数( ,) p Cx X=1,5 X=6, ,10 X=11, ,15 X=16, ,20 .172633e-14 2038.69 .886669e9 .315986e11 .859841e-8

28、 76273.1 .347381e10 .235250e11 .114368e-3 .156701e7 .978210e10 .116152e11 .114955 .195843e8 .199945e11 .341541e10 25.5252 .159475e9 .296669e11 .452548e9 可见而靠近1x处的条件数非常小, 而16x处的条件数最大 . 这与我们的 观测结果是一致的 . 由于条件数太大所以此问题是病态问题. 实验 Y 用 Jacobi法求对称矩阵的特征值及特征向量 一、实验内容 已下列矩阵为例 , 求对称矩阵的全部特征值及特征向量: 1. 5 524 253 ,1

29、10 . 431 A 2. 6 351 529,1 10 . 193 A 3. 9 213 126 ,1 10 . 361 A 二、方法步骤 1. 在A的非主对角线元素中 , 找出按模最大的元素 pq a; 2. 计算平面旋转矩阵 pq U, 其中的sin( )及cos( )由 cot2 2 ppqq pq aa a 计算; 3. 计算 T 1pqpq AUAU, 1pq UU U( 1 U 的初始值取单位阵 ); 4. 如 果 ( 1 ) max ij ij a( 其中 ( 1) ij a为 1 A 的元素 ), 则 停止计 算 , 所求特征值 为: (1) iii a, 特征向量 : T

30、1, , iin i UUx即得第i列; 否则令 1 AA , 1 UU 重复 19 以上各步 . 三、实验结果 讨论的矩阵为 : A = 5 2 4 2 5 3 4 3 1 要求误差为 :err=1e-005 迭代次数 :7 第 1 个特征值为 :9.8057 相应的特征向量为 :0.65053 0.57817 0.49248 第 2 个特征值为 :3.0604 相应的特征向量为 :-0.62538 0.77572 -0.084615 第 3 个特征值为 :-1.8661 相应的特征向量为 :-0.43095 -0.25294 0.8662 讨论的矩阵为 : A = 3 5 -1 5 2 9

31、 -1 9 3 要求误差为 :err=1e-006 迭代次数 :7 第 1 个特征值为 :3.8441 相应的特征向量为 :0.87749 0.052794 -0.47668 第 2 个特征值为 :-8.2319 相应的特征向量为 :-0.3691 0.70896 -0.60095 第 3 个特征值为 :12.3879 相应的特征向量为 :0.30622 0.70327 0.64159 讨论的矩阵为 : A = 2 1 3 1 2 6 3 6 1 要求误差为 :err=1e-009 迭代次数 :8 第 1 个特征值为 :1.1984 相应的特征向量为 :0.88975 -0.4478 -0.088467 第 2 个特征值为 :8.6956 相应的特征向量为 :0.39037 0.64607 0.6559 第 3 个特征值为 :-4.894 相应的特征向量为 :-0.23656 -0.61812 0.74965 以上为程序 :tezhjcobi.m运行结果 ,Matlab 程序见附件 20