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

《数值分析课程设计报告28436.pdf》由会员分享，可在线阅读，更多相关《数值分析课程设计报告28436.pdf（15页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

2、in Matlab environment,including theroy analysis,program and results.Whats 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 Lagranges int

3、erpolation polynomial and Newtons 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.主题词：Romberg 积分，插值方法，Langrange 插值

4、多项式，Newton 插值多项式，拟合多项式。一.问题提出：已知椭圆的周长可以表示成 s=a2022cos1d(01),取 a=1,针对从 0.1 到 0.9（步长 h=0.1）分别求出周长 s;(用 Romberg 积分方法)对于以上数据，求出的插值多项式；对于中数据，试用最小二乘法的思想求作拟合多项式（要求是偶次），并对这些多 项式的优劣进行比较。二问题解决：依据问题出现先后顺序，对三个小问进行讨论。Romberg 算法是在复化梯形求积公式的基础上，应用理查逊外推构造的一种数值积分方法。由复合梯形公式的展开定理，得到如下关系式：T1(h)-I=aah2+a2h4+a3h6+amh2m+其中

5、，I=badxxf)(,T1(h)=Tn.利用 Richardson 外推定理对 T1(h)进行加速，注意这里取 m=1,q=21,有 利用 T0(2h)和 T0(22h)可以得到 实际上 T1(h)就是复化抛物线求积公式，一般的计算公式,)21(1)2()21()2()2()1(21)1(211mkmmkmkmhThThT(m=1,2,k=0,1,2,)即14)2()2(4)(11111mkmkmmmhThThT.由于 Romberg 求积过程是每次把区间缩小一半，所以 Romberg 积分方法也叫做逐次分半加速收敛法。Robermg 求积算法的计算过程如下：取 k=0,h=b-a,求)0(

6、0T=).()(2bfafh 令 1k(k 记区间a,b的二分次数).求梯形值 T0kab2，即按递推公式10212)(221nkknnxfhTT计算)(0kT.求加速值，按公式14)2()2(4)(11111mkmkmmmhThThT逐个求出如表的第 k 行其余 各元素 Tj(k-j)(j=1,2,k).若)0(1)0(kkTT0.0001&ik)|i4 i=i+1;h=h/2;s=0;for n=1:M x=a+h*(2*n-1);s=s+feval(f,p,x);end R(i+1,1)=R(i,1)/2+h*s;M=2*M;for j=1:i R(i+1,j+1)=R(i+1,j)+(

7、R(i+1,j)-R(i,j)/(4j-1);end err=abs(R(i,i)-R(i+1,i+1);end S=R(i+1,i+1)function y=f(p,x)y=sqrt(1+p.2*cos(x).2);在 Matlab 中运行，分别输入 p=0.1,0.9，计算结果列出如表：0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 s 6.2989 6.3456 6.4223 6.5274 6.6592 6.8153 6.9936 7.1919 7.4082 利用的计算结果，求插值多项式。构造 4 阶零矩阵 D 第一列元素R(J,1)存放二分 J 次后的梯形值 利

8、用公式依据 T 表顺序求每行其余元素 R(J,K)，保存在一个特别的下三角矩阵中 当|R(J,J)-R(J+1,J+1)|y=12.6488*power(x,8)-51.4881*power(x,7)+87.9514*power(x,6)-81.7083*power(x,5)+44.6592*power(x,4)-14.8231*power(x,3)+4.3946*power(x,2)-0.2795*x+6.2940 y=Columns 1 through 7 6.2989 6.3456 6.4223 6.5274 6.6592 6.8153 6.9936 Columns 8 through

9、9 7.1920 7.4083 Lagrange 插值多项式形式对称，公式中的每一项与所有的插值点有关。因此，用它计算 f(x)在插值点 x 的近似值时，也有不便之处：假如需要增加一个插值结点，则原来计算的插值多项式 Pn(x)对计算 Pn+1(x)没有用。我们希望，增加新的插值结点时，原先计算的结果对后面的计算过程仍然有用。Newton 插值多项式就具有这种优点。Newton 插值多项式：Nn(x)=a0+a1(x-x0)+an(x-x0)(x-x1)(x-xn-1).其中，ak=fx0,x1,xk,k=0,1,n.一般地，称 111011010,.,.,.,kkkkkkxxxxxfxxxx

10、fxxxf为 f(x)的 k 阶均差。系数 a0,a1,an具有一种不变性。这是说，当计算得到 Nn(x)后，若要增加一个结点，则只要从 Nn(x)增加一项 an+1(x-x0)(x-x1)(x-xn-1)(x-xn)，便立刻得到 Nn+1(x).Newton 插值多项式还可写作：Nn(x)=(an(x-xn-1)+an-1(x-xn-2)+)+a0*均差表如下：xk fxk f,f,f,f,x0 fx0 x1 fx1 fx0,x1 x2 fx2 fx1,x2 fx0,x1,x2 x3 fx3 fx2,x3 fx1,x2,x3 fx0,x1,x2,x3 x4 fx4 fx3,x4 fx2,x3

11、,x4 fx1,x2,x3,x4 fx0,x1,x2,x3,x4 程序框图：构造n阶零矩阵D 把y的值存入 D的第一列 依次计算第 k 阶均差，保留在 D 的第(k+1)列 编写 Matlab 程序 newpoly.m：function C,D=newpoly(x,y)n=length(x);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

12、(m)=C(m)+D(k,k);end 运行结果：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.7830 1.1100 -0.3500 0.1250 0.8333 2.2222 7.1919 1

13、.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.3946 Columns 8 through 9 根据(*)式利用差商即 D对角线上值元素求得Newton插值

14、多项式 -0.2795 6.2940 将结果整理，得到 Newton 插值多项式为：y=12.6488x8-51.4881x7+87.9514x6-81.7083x5+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.8153 6.9936 7.1919 7.4082;c

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

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

17、小误差平方和，表示为：为求解 E 的最小值，每个偏导数必须为零（即MicEi,.,2,1,0/），这样可得到如下方程组：交换上式中的求和顺序，可得一个 M*M 线性方程组，未知数是系数cj，称为正规方程组：将上面的正规方程组表示成矩阵形式可减少不必要的计算量。方程组*的系数矩阵是一个对称矩阵，并且是正定的，可以唯一地确定系数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

18、(x2)f3(xN)fM(x1)fM(x2)fMxN)乘积 FF 是一个 M*M 矩阵，方程组*转化为求解线性方程组：YFFCF（求解系数矩阵 C）程序框图：编写 Matlab 程序 lspoly.m:function c=lspoly(x,y,M)n=length(x);B=zeros(1:M+1);F=zeros(n,M+1);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 6.8153 6.9936 7.1919 7.408

19、2;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.4769 Column 9 12.6460 整理得到偶次拟合多项式：2 阶拟合多

20、项式：y=6.2624+0.1864x+1.2147x2 构造 M+1 维零向量 B，N 行 M+1 列的零矩阵 F 将 xk-1依次存放到 F的第 k 列中 求解正规方程组 FFC=FY 4 阶拟合多项式：y=6.2844-0.02x+1.6838x2-0.2704x3-0.0382x4 6 阶拟合多项式：y=6.2829+0.0069x+1.5243x2+0.1506x3-0.5697x4+0.3013x5-0.0556x6 8 阶拟合多项式：y=6.2940-0.2795x+4.3940 x2-14.8201x3+44.6498x4-81.6911x5+87.9326x6-51.4769

21、x7+12.6460 x8 为检验每个拟合多项式的效果，编写程序 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);将结果列出如表：实际值 6.2989 6.3456 6.4223 6.5274 6.6592 6.8153 6.9936 7

22、.1919 7.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 阶拟合多项式的拟合结果已与实

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

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

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

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

27、比较好，因为高阶多项式的拟合会引起数值的不稳定，使用最小二乘法多项式拟合非线性数据的方法固然具有吸引力，但如果数据不具有多项式特性，则求出的曲线可能产生大的振荡。这种现象在高阶多项式情况下更容易发生，由于这个原因，一般很少使用超过 6 阶的多项式。正因为这样，所以才会选 6 阶拟合多项式。参考文献：John H.Mathews,Kurtis D.Fink 著，陈渝等译。Numerical Methods Using MATLAB(Third Edition)北京：电子工业出版社.张池平主编.计算方法.哈尔滨：哈尔滨工业大学出版社.王兵团，桂文豪著.数学实验基础.北京：北方交通大学出版社.杨万利

28、等著.数值分析教程.北京：国防工业出版社.康晓东主编.数值算法与非数值算法.北京：电子工业出版社 徐士良编著.数值分析与算法.北京：机械工业出版社.林成森编著.数值计算方法.北京：科学出版社.张池平，施云慧著.计算方法.北京：科学出版社.西安郵電學院 数值分析课程设计报告书 系 部 名 称：应 用 数 理 系 学 生 姓 名：专 业 名 称：信息与计算科学 班 级：信息 070 x 时间：2 实验 X 19 次方系数对方程根的影响 一、问题的提出 考虑代数多项式201()(1)(2)(20)()kP xxxxxk,显然它的全部根为 1,2,20;如果19次方的系数有扰动,则对方程的解有何影响.

29、二、实验内容 直接利用 Matlab 中的 roots()和 poly()函数,选择不同的扰动做数值计算.主要步骤:ve=zeros(1,21);ve(2)=ess;roots(poly(1:20)+ve).Matlab 程序:见附件 bingtaiwenti.m 三、实验结果 1.取扰动分别为36910,10,10所得结果如下:310 610 910 27.0817+5.03812i 21.3025+1.56717i 19.952+0i 19.2293+0i 19.5337+9.1664i 18.5028+3.6004i 17.6573+0.692896i 13.8235+7.77167i

30、15.1651+3.76125i 15.4524+0.875524i 10.7211+5.4609i 12.4866+2.88278i 13.3527+0.486992i 8.91282+3.47317i 10.5225+1.71958i 11.8578+0i 11.0427+0i 7.69268+1.89884i 9.04487+0.59455i 9.9916+0i 9.00201+0i 6.75761+0.654714i 7.94891+0i 7.99952+0i 7.00247+0i 7.00009+0i 5.95208+0i 5.99995+0i 5.99999+0i 5.00061+

31、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 2.绘制相应的图形,即将以上结果可视化.四、实验结果分析 观测现象:由结果可以观察到误差是不可避免的,并且扰动的减小不能使得所有根除的误差减小,还可以观察到16x 处的误差总是最大,而靠近1x 处的误差较小.误差分析:原问题数学上描述为:19()0P xx的求解问题,也可等价为:则要考虑的问题是对方程的某个解x的影响,我们不妨将x看成的函数,从而可以求得此问题的条件数为:可以求得20 个根处的条件数为:19 次方系数变化对根的影响问题条件数(,)pC

32、x X=1,5 X=6,10 X=11,15 X=16,20.172633e-14 2038.69.886669e9.315986e11.859841e-8 76273.1.347381e10.235250e11.114368e-3.156701e7.978210e10.116152e11.114955.195843e8.199945e11.341541e10 25.5252.159475e9.296669e11.452548e9 可见而靠近1x 处的条件数非常小,而16x 处的条件数最大.这与我们的观测结果是一致的.由于条件数太大所以此问题是病态问题.实验Y 用Jacobi法求对称矩阵的特征

33、值及特征向量 一、实验内容 已下列矩阵为例,求对称矩阵的全部特征值及特征向量:1.5524253,1 10.431A 2.6351529,1 10.193A 3.9213126,1 10.361A 二、方法步骤 1.在A的非主对角线元素中,找出按模最大的元素pqa;2.计算平面旋转矩阵pqU,其中的sin()及cos()由cot22ppqqpqaaa计算;3.计算T1pqpqAUAU,1pqUU U(1U的初始值取单位阵);4.如果(1)maxijija(其中(1)ija为1A的元素),则停止计算,所求特征值为:(1)iiia,特征向量:T1,iin iUUx即得第i列;否则令1AA,1UU重

34、复以上各步.三、实验结果 讨论的矩阵为: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 -1 9 3 要求误差为:err=1e-006 迭代次数:7 第 1 个特征值为:3.8441 相应的特征向量为:0.

35、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 程序见附件