《数值分析课程设计报告》由会员分享，可在线阅读，更多相关《数值分析课程设计报告（20页珍藏版）》请在金锄头文库上搜索。
1、1数 值 分 析 课 程 设 计 报 告摘要（中）：本文建立在数值分析的理论基础上，能够在 Matlab 环境中运行，给出了理论分析、程序清单以及计算结果。更重要的是， 还有详细的对算法的框图说明。首先运用 Romberg 积分方法对给出定积分进行积分，然后对得到的结果用插值方法，分别求出 Lagrange 插值多项式和 Newton 插值多项式，再运用最小二乘法的思想求出拟合多项式，最后对这些不同类型多项式进行比较，找出它们各自的优劣。Summary(Enlish):This essay is based on numerical analysis ,and could be operate
2、d 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 Lag
3、ranges interpolation 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 积分，插值
4、方法， Langrange 插值多项式，Newton 插值多项式，拟合多项式。2一. 问题提出：已知椭圆的周长可以表示成 s=a (00.0001&iy=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.2940y =Columns 1 through 7 6.2989 6.3456 6.4223 6.5274 6.6592 6.8153 6.9936Colum
5、ns 8 through 9 7.1920 7.4083Lagrange 插值多项式形式对称，公式中的每一项与所有的插值点有关。因此，用它计算 f(x)在插值点 x 的近似值时，也有不便之处：假如需要增加一个插值结点，则原来计算的插值多项式 Pn(x)对计算 Pn+1(x)没有用。我们希望，增加新的插值结点时，原先计算的结果对后面的计算过程仍然有用。Newton 插值多项式就具有这种优点。 Newton 插值多项式：N n(x)=a0+a1(x-x0)+an(x-x0)(x-x1)(x-xn-1).其中，a k=fx0,x1,xk,k=0,1,n.一般地，称为 f(x)的 k 阶均差。1101
6、010 ,.,.,. kkkk xxfxfxf系数 a0,a1,an 具有一种不变性。这是说，当计算得到 Nn(x)后，若要增加一个结7点，则只要从 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 fx0x1 fx1 fx0,x1x2 fx2 fx1,x2 fx0,x1,x2x3 fx3 fx2,x3 fx1,x2,x3 fx0,x
7、1,x2,x3x4 fx4 fx3,x4 fx2,x3,x4 fx1,x2,x3,x4 fx0,x1,x2,x3,x4程序框图：编写 Matlab 程序 newpoly.m：function C,D=newpoly(x,y)n=length(x);D=zeros(n,n);D(:,1)=y;for j=2:nfor k=j:nD(k,j)=(D(k,j-1)-D(k-1,j-1)/(x(k)-x(k-j+1);endendDC=D(n,n);for k=(n-1):-1:1C=conv(C,poly(x(k);m=length(C);C(m)=C(m)+D(k,k);end构造 n 阶零矩阵D
8、把 y 的值存入 D 的第一列依次计算第 k 阶均差，保留在 D 的第(k+1) 列中根据(*)式利用差商即 D 对角线上值元素求得 Newton 插值多项式8运行结果：D =Columns 1 through 7 6.2989 0 0 0 0 0 06.3456 0.4670 0 0 0 0 06.4223 0.7670 1.5000 0 0 0 06.5274 1.0510 1.4200 -0.2667 0 0 06.6592 1.3180 1.3350 -0.2833 -0.0417 0 06.8153 1.5610 1.2150 -0.4000 -0.2917 -0.5000 06.9
9、936 1.7830 1.1100 -0.3500 0.1250 0.8333 2.22227.1919 1.9830 1.0000 -0.3667 -0.0417 -0.3333 -1.94447.4082 2.1630 0.9000 -0.3333 0.0833 0.2500 0.9722Columns 8 through 9 0 00 00 00 00 00 00 0-5.9524 04.1667 12.6488ans =Columns 1 through 7 12.6488 -51.4881 87.9514 -81.7083 44.6592 -14.8231 4.3946Columns
10、 8 through 9 -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.mfunction newpolyplotx=0.1:0.1:0.9;y=6.2989 6.3456 6.4223 6.5274 6.6592 6.8153 6.9936 7.1919 7.4
11、082;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 插值那样的均差表，在插值点变化时，有些数据可以保留。 拟合多项式曲线拟合问题与函数插值问题不同，它不要求曲线通过所有
12、已知点，只要求得到的近似函数能反映数据的基本关系。因此，曲线拟合的过程比插值过程得到的结果更能反映客观实际。在某种意义上，曲线拟合更具有实用价值，因为实际问题中所提供的数据往往很多，如果用插值法势必要得到次数很高的插值多项式，导致计算上的很多麻烦。一般总是希望各观察数据与拟合曲线的偏差的平方和最小，这样就能使拟合曲线更接近于真实函数。这个原理就称为最小二乘原理。用最小二乘原理作为衡量“曲线拟合优劣”的准则称为曲线的最小二乘法。最小二乘法的思想是：设有 N 个数据点 （x k,yk），并给定 M 个线性独立函数fj(x)。为求 M 个稀疏，使用由线性组合形成的函数 f(x)，表示为：jjxfcx
13、f1求解最小误差平方和，表示为： Nk kMjjMkk yxfcyxfcE121221,.为求解 E 的最小值，每个偏导数必须为零（即 ） ，这样MiEi,.1,0/10可得到如下方程组： MixfyxfcNk kikMjj ,.21,01 交换上式中的求和顺序，可得一个 M*M 线性方程组，未知数是系数c j，称为正规方程组： Mj NkkijNkkjj iyxfcxf1 1,.21, *将上面的正规方程组表示成矩阵形式可减少不必要的计算量。方程组 的系数矩阵是一个对称矩阵，并且是正定的，可以唯一地确定系数c j，*根据这些特点，构造如下矩阵：f1(x1) f2(x1) fM(x1)f1(x
14、2) 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)乘积 FF 是一个 M*M 矩阵，方程组 转化为求解线性方程组：*（求解系数矩阵 C）YFC程序框图：编写 Matlab 程序 lspoly.m:function c=lspoly(x,y,M)n=length(x);B=zeros(1:M+1);构造 M+1 维零向量 B，N 行 M+1 列的零矩阵 F将 xk-1 依次存放到 F 的第 k 列中求解正规方程组 FFC=FY11F=zeros(n,M+1);for k=1:M+1F(:,k)=x.(k-1);endA=F*F;B=F*y;c=(AB);运行结果为： x=0.1:0.1:0.9; y=6.2989 6.3456 6.4223 6.5274 6.6
