数值分析插值上机报告含代码Word格式文档下载.docx

1、 多项式插值一般有两种方法，拉格朗日插值和牛顿插值。个人认为，牛顿插值便于计算机实现，所以采用牛顿插值法来解决第一问。但是牛顿插值公式并不十分直观，最好能将牛顿插值多项式化简。为了化简多项式，我建立了一个多项式类，在进行多项加法和乘法时可以自动化简多项式，在逐步建立牛顿插值多项式时，该多项式会自动化为最简。题目中还要 求在同一张图中画出原多项式和插值多项式，可以使用MFC进行画图，通过菜单选项来选择等分点的数量，然后程序会自动算出插值多项式并输出在屏幕上，接着会输出图形。 第二问： 第二问要求求出最大相对误差，本题较为特殊，其误差表达式恰好为一个多项式，所以求误差多项式后，逐点搜索该函数的绝对

2、值的最大值。 第三问： 该问要求分析前两问所得结果，同时要求寻找更好的方法来插值。本人使用了三次样条插值（边界条件为一阶导数相等）来求解该问题。求三次样条函数的关键过程其实就是求解三对角方程，可以使用追赶法来求解该方程。为了体现效果，也可以用MFC来画图，同时输出各段插值函数以及最大相对误差。 求解方法： 问题中用到的两个插值方法为牛顿插值和三次样条插值，解法与教科书中的一模一样，在此不再赘述。 程序设计： 1.牛顿插值 程序设计分为三步走： 设计多项式类 在附录一的poly.h中 设计牛顿插值函数 在附录一的 Newton_interpolation.h中 设计MFC程序 关键代码代码在附录

3、一的 MFC中 多项式类的声明如下： typedef int degree; typedef double coefficient; /term represents polynomial term consisting of its power and coefficient typedef std:pair term; class polynomial private : /liner list std:list poly; public: /default constructor polynomial（）; /copy constructor polynomial（const polyn

4、omial / construct a polynomial consisting of only one term coefficient*Xdegree 从声明中可以看到，该类使用STL中的list类型来存储多项式，还实现了多项式的加法和乘法，因为这两个运算在插值过程中大量使用。类中最后一个操作是重载了（）运算符，该操作是计算多项式在x上的值，其实是为了计算最大相对误差和画图之用。 牛顿插值函数的过程大致为： 求出差商表 计算出插值公式 MFC设计过程中关键步骤有二，一为建立菜单供选择等分点数数目之用，在菜单处理函数中设置等分点个数并调用视图类的操作InValidate（）使当前画面无效，

5、系统再自动调用OnDraw（）重绘图形。其二为设计OnDraw（）函数，该函数为视图类的一个操作用于绘图，用户可以重载，在该函数中，输出插值函数，画出坐标系以及图形。 2.三次样条插值 三次样条插值函数基于以上的多项式类和线性方程组中的追赶法。整个实现较为简单，在附录一中的Sample.h中。为了考察插值效果，将样条插值的结果显示在MFC中，整个架构是也是通过菜单来选择等分点的数目，再显示图形及具体函数。与牛顿插值一起整合在MFC中。 小程序的截图如下： 该图是牛顿插值法N=10时的情况，可通过菜单选择N的值 该图是三次样条插值N=10的情况，可通过菜单选择N的值 结果分析： 牛顿插值法的图如

6、下： N=2 N=4 N=6 N=8 N=10 可以看到以下几点 1. 随着N的增大，插值多项式的波动愈发剧烈，最后，除了插值点外，函数形状与原函数形状完全偏离。此时发生了龙格（Runge）现象。 2. 如果仔细观察最大相对误差（图中字太小，最好运行程序看一下）的话，发现误差不断增大，而且增速是变快的。 3. 还可以看到，插值多项式只有偶次项，因为原函数函数是偶函数。 4. 越偏离原点，插值函数与原函数的偏差越大。这一点是多项式函数与原函数的性质造成的，原函数随着X的绝对值趋近无穷大而趋向于零，而多项式随着X的绝对值趋近无穷大而趋向于无穷大，所以两者的偏差会随着|X|的变大而变大。这一点也反映

7、在最大相对误差往往出现的靠近与?1的地方。 由此看出，多项式插值时不宜选取次数太高的多项式，为了避免Runge现象，有效途径之一就是分段插值，其中效果比较好的是三次样条插值。 三次样条插值的图如下： N=10 从上方五张图可以看到随着N的增大，样条插值多项式越来越接近原函数，到N=10的时候，两条曲线几乎难以分辨。如果运行程序仔细看的话，最大相对误差也是逐渐趋向于零，实际上当N=10的时候，最大相对误差为0.033左右。但是插值公式也越来越复杂，不易使用。 vi=0.5/（2-0.5*vi）; ui=（di-ui-1*0.5）/（2-0.5*vi-1）; for（int i=n-1;ii-）

8、mi=ui-vi*mi+1; m=dfp（low）;mn=dfp（up）; for（int i=1;ii+） polynomial temp（0,fp（low+h*i-h）/h/h/h）; temp*=（polynomial（0,-（low+h*i）+polynomial（1,1）; temp*=（polynomial（1,2）+polynomial（0,h-2*（low+h*i-h）; Si+=temp; temp=polynomial（0,fp（low+h*i）/h/h/h）; temp*=（polynomial（0,-（low+h*i-h）+polynomial（1,1）; temp*=

9、（polynomial（1,-2）+polynomial（0,h+2*（low+h*i）; temp=polynomial（0,mi-1/h/h）; temp*=（polynomial（1,-1）+polynomial（0,low+h*i）; temp*=（polynomial（1,1）+polynomial（0,-（low+h*i-h）; temp=polynomial（0,mi/h/h）; temp*=（polynomial（1,1）+polynomial（0,-（low+h*i）; MFC void CMFCView:OnDraw（CDC* pDC） CMFCDoc* pDoc = Ge

10、tDocument（）; ASSERT_VALID（pDoc）; if （!pDoc） return; if（Newton） polynomial p1; p1=Newton_interpolation（f2,N,-1,1）; CString str; std:string ss;ostringstream strm; strm TextOutW（0,0,str）; pDC-TextOut（0,25,str2）; pDC-SetTextColor（RGB（255,0,0）; str3.Format（L最大相对误差在 X=%f 或 %f 时 取得，RelativeError=%f ,p,-p,m

11、ax）; pDC-TextOut（0,50,str3）;/p1; pDC-SetTextColor（RGB（0,0,0）; CRect rect; GetClientRect（rect）; pDC-SetViewportOrg（rect.Width（）/2,rect.Height（）*3/4）; CRect rect1（420,110,-420,-410）; pDC-Rectangle（rect1）; pDC-MoveTo（-420,0）; pDC-LineTo（420,0）; pDC-MoveTo（0,110）; pDC-LineTo（0,-410）; for（double d=-1;dd+

12、=0.1） int posx=d*400,posy=0; pDC-MoveTo（posx,posy）; pDC-LineTo（posx,posy+4）; str.Format（L,d）; pDC-TextOutW（posx,posy+6,str）; for（double d=-0.5;dd+=0.1） if（abs（d）0.01） int posx=0,posy=-200*d; pDC-LineTo（posx+4,posy）; pDC-TextOutW（posx+6,posy,str）; CPen redpen,bluepen,greenpen; CBrush greenbrush; redp

13、en.CreatePen（PS_SOLID,2,RGB（255,0,0）; bluepen.CreatePen（PS_SOLID,2,RGB（0,0,255）; greenpen.CreatePen（PS_SOLID,1,RGB（0,255,0）; greenbrush.CreateSolidBrush（RGB（0,255,0）; pDC-SelectObject（redpen）; double x,y; for（int i=0;ii+） xi=-1+0.02*i; yi=-200/（1+25*xi*xi）; xi*=400; pDC-MoveTo（x,y）;ii+） pDC-LineTo（x

14、i,yi）; pDC-MoveTo（xi,yi）; pDC-SelectObject（bluepen）; yi=-200*p1（xi）; pDC-SelectObject（greenpen）; pDC-SelectObject（greenbrush）; for（double h=2.0/N,d=-1;d1.05;d+=h） int posx=400*d,posy=-200*f2（d）; pDC-Ellipse（CRect（posx-3,posy-3,posx+3,posy+3）; else CString str; double h=2.0/N; str.Format（L当N=%d时 三次样条

15、的各段插值函数为,N）; pDC-TextOutW（0,0,str）; vector S（N+1）; Sample（f2,df,N,-1,1,S）;ii+） std: strm TextOut（0,20*i,str2）; CRect rect; pDC-SetViewportOrg（rect.Width（）/2+250,rect.Height（）*3/4）; for（double d=-0.5; CPen redpen,bluepen,greenpen;/si;/polynomial double maxpos=-1,maxerr=0,temp; temp=Si*N/101+1（xi）; if

16、（abs（temp-f2（xi）/f2（xi）maxerr） maxpos=xi;maxerr=abs（temp-f2（xi）/f2（xi）; yi=-200*temp; pDC-SelectObject（greenpen）; pDC-SetViewportOrg（0,20*N+20）; str.Format（L最大相对误差在 ）; str.AppendFormat（L或 %4.2f时取得,为 %f,maxpos,-maxpos,maxerr）; Matrix.h #pragma once #include cmathiomanip template class Matrix /iomanip

17、/cmath/iostream private: T *mat; int row; int column; / used in destructor and assignment operator void destroy（）; /constructor destructor and assignment /default constructor Matrix（）:mat（NULL）,row（0）,column（0） ; /construct a n*n identity matrix explicit Matrix（int n）; /construct a m*n zero matrix M

18、atrix（int m,int n, T def=0）; / copy constructor Matrix（const Matrix / assignment operator Matrix operator=（const Matrix /destructor Matrix（） destroy（）; /miscellaneous procedures / return the number of row int Row（） return row; / return the number of column int Column（） return column; / return the re

19、ference of Mmn.This kind of grammar is like matlab. T operator（） （int m,int n）; / overload the ostream template friend std:ostream operator（std:ostream os,Matrix h）; ;/type template void Matrix :destroy（） if（mat!=NULL） for（int i=0;i /t Matrix Matrix（int n）:mat（NULL）,row（0）,column（0） if（n=0） std:cerr

20、 std:endl; else mat=new t *n; for（int i=in;i+） mati= tn; matii=1; /std:i+） for（int j=jj+） if（i!=j） matij=/n;Matrix（int m,int n, T def）:mat（NULL）,row（0）,column（0） if（m0 n0） row=m; column=n; mat=new T *m;m;i+） mati= j=def;/m;Matrix（const Matrix M）:mat（NULL）,row（0）,column（0） if（M.mat!=NULL） mat=new T *M.row;m.row;i+） mati= tm.column; for（int j=M.column;M.matij; row=M.row;/m.row;operator=（const Matrix M） if（this!=M） destroy（）; if（M.mat! else mat=NULL; column=0; retur