数值分析插值上机报告含代码.docx

数值分析插值上机报告含代码.docx

《数值分析插值上机报告含代码.docx》由会员分享，可在线阅读，更多相关《数值分析插值上机报告含代码.docx（8页珍藏版）》请在冰豆网上搜索。

数值分析插值上机报告含代码

数值分析插值上机报告含代码

许多实际问题的内在规律都可以用函数y?

f（x）进行确切描述，但是往往该函数的确切表达式是不可知的，容易得到的是函数的一些采样点，我们有时希望观察函数的大致形状或者希望知道函数的解析式，此时可以借助插值的方法的来求得函数的大致表达式。

但是插值法毕竟是一种近似方法，有一定的误差，所以我们需要分析什么时候插值法是不适合的，以及何种插值法比较适合某个具体问题。

借助书中的实验问题，本文考察了多项式插值的缺陷，画出了在各个情况下插值曲线和原函数曲线的图像，并且计算了最大的相对误差，发现随着多项式的次数的增大，最大相对误差不断增大。

还实现了三次样条插值法，看到三次样条插值的效果非常好，并且随着等分点的数量的增大，最大相对误差不断减小。

问题重述：

对于函数

f（x）?

1x?

[?

1,1]

1?

（5x）22。

n取n+1个基点，xi?

?

1?

ih（i?

0,1,2,?

n），其中h?

（1）对n?

2,4,6,8,10分别作n次插值多项式Pn（x），并在同一坐标系分别画出f（x）与Pn（x）。

（2）在非节点处分别计算f（x）与Pn（x）的最大相对误差：

?

1?

x?

1x?

ximax?

f（x）?

Pn（x）n?

2,4,6,8,10

f（x）（3）根据f（x）与Pn（x）的图形及最大相对误差进行比较分析，试寻找插值效果较好的改进方法。

问题分析：

第一问：

多项式插值一般有两种方法，拉格朗日插值和牛顿插值。

个人认为，牛顿插值便于计算机实现，所以采用牛顿插值法来解决第一问。

但是牛顿插值公式并不十分直观，最好能将牛顿插值多项式化简。

为了化简多项式，我建立了一个多项式类，在进行多项加法和乘法时可以自动化简多项式，在逐步建立牛顿插值多项式时，该多项式会自动化为最简。

题目中还要

求在同一张图中画出原多项式和插值多项式，可以使用MFC进行画图，通过菜单选项来选择等分点的数量，然后程序会自动算出插值多项式并输出在屏幕上，接着会输出图形。

第二问：

第二问要求求出最大相对误差，本题较为特殊，其误差表达式恰好为一个多项式，所以求误差多项式后，逐点搜索该函数的绝对值的最大值。

第三问：

该问要求分析前两问所得结果，同时要求寻找更好的方法来插值。

本人使用了三次样条插值（边界条件为一阶导数相等）来求解该问题。

求三次样条函数的关键过程其实就是求解三对角方程，可以使用追赶法来求解该方程。

为了体现效果，也可以用MFC来画图，同时输出各段插值函数以及最大相对误差。

求解方法：

问题中用到的两个插值方法为牛顿插值和三次样条插值，解法与教科书中的一模一样，在此不再赘述。

程序设计：

1.牛顿插值

程序设计分为三步走：

设计多项式类在附录一的poly.h中设计牛顿插值函数在附录一的Newton_interpolation.h中设计MFC程序关键代码代码在附录一的MFC中

多项式类的声明如下：

typedefintdegree;typedefdoublecoefficient;//termrepresentspolynomialtermconsistingofitspowerandcoefficienttypedefstd:

:

pairterm;classpolynomial{private:

//linerliststd:

:

listpoly;public:

//defaultconstructorpolynomial（）{};//copyconstructorpolynomial（constpolynomial//constructapolynomialconsistingofonlyonetermcoefficient*X^degree

从声明中可以看到，该类使用STL中的list类型来存储多项式，还实现了多项式的加法和乘法，因为这两个运算在插值过程中大量使用。

类中最后一个操作是重载了（）运算符，该操作是计算多项式在x上的值，其实是为了计算最大相对误差和画图之用。

牛顿插值函数的过程大致为：

求出差商表

计算出插值公式

MFC设计过程中关键步骤有二，一为建立菜单供选择等分点数数目之用，在菜单处理函数中设置等分点个数并调用视图类的操作InValidate（）使当前画面无效，系统再自动调用OnDraw（）重绘图形。

其二为设计OnDraw（）函数，该函数为视图类的一个操作用于绘图，用户可以重载，在该函数中，输出插值函数，画出坐标系以及图形。

2.三次样条插值

三次样条插值函数基于以上的多项式类和线性方程组中的追赶法。

整个实现较为简单，在附录一中的Sample.h中。

为了考察插值效果，将样条插值的结果显示在MFC中，整个架构是也是通过菜单来选择等分点的数目，再显示图形及具体函数。

与牛顿插值一起整合在MFC中。

小程序的截图如下：

该图是牛顿插值法N=10时的情况，可通过菜单选择N的值

该图是三次样条插值N=10的情况，可通过菜单选择N的值

结果分析：

牛顿插值法的图如下：

N=2

N=4

N=6

N=8

N=10可以看到以下几点

1.随着N的增大，插值多项式的波动愈发剧烈，最后，除了插值点外，函数形状与原函数形状完全偏离。

此时发生了龙格（Runge）现象。

2.如果仔细观察最大相对误差（图中字太小，最好运行程序看一下）的话，发现误差不断增大，而且增速是变快的。

3.还可以看到，插值多项式只有偶次项，因为原函数函数是偶函数。

4.越偏离原点，插值函数与原函数的偏差越大。

这一点是多项式函数与原函数的性质造成的，原函数随着X的绝对值趋近无穷大而趋向于零，而多项式随着X的绝对值趋近无穷大而趋向于无穷大，所以两者的偏差会随着|X|的变大而变大。

这一点也反映在最大相对误差往往出现的靠近与?

1的地方。

由此看出，多项式插值时不宜选取次数太高的多项式，为了避免Runge现象，有效途径之一就是分段插值，其中效果比较好的是三次样条插值。

三次样条插值的图如下：

N=2

N=4

N=6

N=8

N=10

从上方五张图可以看到随着N的增大，样条插值多项式越来越接近原函数，到N=10的时候，两条曲线几乎难以分辨。

如果运行程序仔细看的话，最大相对误差也是逐渐趋向于零，实际上当N=10的时候，最大相对误差为0.033左右。

但是插值公式也越来越复杂，不易使用。

}

{v[i]=0.5/（2-0.5*v[i]）;u[i]=（d[i]-u[i-1]*0.5）/（2-0.5*v[i-1]）;}

for（inti=n-1;ii--）m[i]=u[i]-v[i]*m[i+1];m=dfp（low）;m[n]=dfp（up）;for（inti=1;ii++）{polynomialtemp（0,fp（low+h*i-h）/h/h/h）;temp*=（polynomial（0,-（low+h*i））+polynomial（1,1））;temp*=（polynomial（0,-（low+h*i））+polynomial（1,1））;temp*=（polynomial（1,2）+polynomial（0,h-2*（low+h*i-h）））;S[i]+=temp;temp=polynomial（0,fp（low+h*i）/h/h/h）;temp*=（polynomial（0,-（low+h*i-h））+polynomial（1,1））;temp*=（polynomial（0,-（low+h*i-h））+polynomial（1,1））;temp*=（polynomial（1,-2）+polynomial（0,h+2*（low+h*i）））;S[i]+=temp;temp=polynomial（0,m[i-1]/h/h）;temp*=（polynomial（1,-1）+polynomial（0,low+h*i））;temp*=（polynomial（1,-1）+polynomial（0,low+h*i））;temp*=（polynomial（1,1）+polynomial（0,-（low+h*i-h）））;S[i]+=temp;temp=polynomial（0,m[i]/h/h）;temp*=（polynomial（1,1）+polynomial（0,-（low+h*i-h）））;temp*=（polynomial（1,1）+polynomial（0,-（low+h*i-h）））;temp*=（polynomial（1,1）+polynomial（0,-（low+h*i）））;S[i]+=temp;}

MFC

voidCMFCView:

:

OnDraw（CDC*pDC）{CMFCDoc*pDoc=GetDocument（）;ASSERT_VALID（pDoc）;if（!

pDoc）return;if（Newton）{polynomialp1;p1=Newton_interpolation（f2,N,-1,1）;CStringstr;std:

:

stringss;std:

:

ostringstreamstrm;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,max）;pDC-TextOut（0,50,str3）;

pDC-SetTextColor（RGB（0,0,0））;CRectrect;

GetClientRect（rect）;

pDC-SetViewportOrg（rect.Width（）/2,rect.Height（）*3/4）;CRectrect1（420,110,-420,-410）;pDC-Rectangle（rect1）;pDC-MoveTo（-420,0）;pDC-LineTo（420,0）;pDC-MoveTo（0,110）;pDC-LineTo（0,-410）;

for（doubled=-1;dd+=0.1）{intposx=d*400,posy=0;pDC-MoveTo（posx,posy）;pDC-LineTo（posx,posy+4）;str.Format（L\,d）;pDC-TextOutW（posx,posy+6,str）;}

for（doubled=-0.5;dd+=0.1）{if（abs（d）0.01）{intposx=0,posy=-200*d;pDC-MoveTo（posx,posy）;pDC-LineTo（posx+4,posy）;str.Format（L\,d）;pDC-TextOutW（posx+6,posy,str）;}}

CPenredpen,bluepen,greenpen;CBrushgreenbrush;

redpen.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）;doublex,y;for（inti=0;ii++）{x[i]=-1+0.02*i;y[i]=-200/（1+25*x[i]*x[i]）;x[i]*=400;}

pDC-MoveTo（x,y）;for（inti=0;ii++）{pDC-LineTo（x[i],y[i]）;pDC-MoveTo（x[i],y[i]）;}

pDC-SelectObject（bluepen）;for（inti=0;ii++）{x[i]=-1+0.02*i;y[i]=-200*p1（x[i]）;x[i]*=400;}

pDC-MoveTo（x,y）;for（inti=0;ii++）{pDC-LineTo（x[i],y[i]）;

pDC-MoveTo（x[i],y[i]）;}pDC-SelectObject（greenpen）;pDC-SelectObject（greenbrush）;for（doubleh=2.0/N,d=-1;d1.05;d+=h）{intposx=400*d,posy=-200*f2（d）;pDC-Ellipse（CRect（posx-3,posy-3,posx+3,posy+3））;}}else{CStringstr;doubleh=2.0/N;str.Format（L\当N=%d时三次样条的各段插值函数为\,N）;pDC-TextOutW（0,0,str）;vector

S（N+1）;Sample（f2,df,N,-1,1,S）;for（inti=1;ii++）{std:

:

stringss;std:

:

ostringstreamstrm;strmTextOut（0,20*i,str2）;}CRectrect;GetClientRect（rect）;pDC-SetViewportOrg（rect.Width（）/2+250,rect.Height（）*3/4）;CRectrect1（420,110,-420,-410）;pDC-Rectangle（rect1）;pDC-MoveTo（-420,0）;pDC-LineTo（420,0）;pDC-MoveTo（0,110）;pDC-LineTo（0,-410）;for（doubled=-1;dd+=0.1）{intposx=d*400,posy=0;pDC-MoveTo（posx,posy）;pDC-LineTo（posx,posy+4）;str.Format（L\,d）;pDC-TextOutW（posx,posy+6,str）;}for（doubled=-0.5;dd+=0.1）{if（abs（d）0.01）{intposx=0,posy=-200*d;pDC-MoveTo（posx,posy）;pDC-LineTo（posx+4,posy）;str.Format（L\,d）;pDC-TextOutW（posx+6,posy,str）;}}CPenredpen,bluepen,greenpen;CBrushgreenbrush;redpen.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）;doublex,y;for（inti=0;ii++）{x[i]=-1+0.02*i;y[i]=-200/（1+25*x[i]*x[i]）;x[i]*=400;}

pDC-MoveTo（x,y）;for（inti=0;ii++）{pDC-LineTo（x[i],y[i]）;pDC-MoveTo（x[i],y[i]）;}

pDC-SelectObject（bluepen）;doublemaxpos=-1,maxerr=0,temp;for（inti=0;ii++）{x[i]=-1+0.02*i;temp=S[i*N/101+1]（x[i]）;if（abs（temp-f2（x[i]））/f2（x[i]）maxerr）{maxpos=x[i];maxerr=abs（temp-f2（x[i]））/f2（x[i]）;}y[i]=-200*temp;x[i]*=400;}

pDC-MoveTo（x,y）;for（inti=0;ii++）{pDC-LineTo（x[i],y[i]）;pDC-MoveTo（x[i],y[i]）;}

pDC-SelectObject（greenpen）;pDC-SelectObject（greenbrush）;

for（doubleh=2.0/N,d=-1;d1.05;d+=h）{intposx=400*d,posy=-200*f2（d）;pDC-Ellipse（CRect（posx-3,posy-3,posx+3,posy+3））;}

pDC-SetViewportOrg（0,20*N+20）;pDC-SetTextColor（RGB（255,0,0））;str.Format（L\最大相对误差在\）;

str.AppendFormat（L\或%4.2f时取得,为%f\,maxpos,-maxpos,maxerr）;pDC-TextOutW（0,0,str）;

Matrix.h

#pragmaonce

#include#include#includetemplateclassMatrix

{

private:

T**mat;introw;intcolumn;//usedindestructorandassignmentoperatorvoiddestroy（）;public:

//constructordestructorandassignment//defaultconstructorMatrix（）:

mat（NULL）,row（0）,column（0）{};//constructan*nidentitymatrixexplicitMatrix（intn）;//constructam*nzeromatrixMatrix（intm,intn,Tdef=0）;//copyconstructorMatrix（constMatrix//assignmentoperatorMatrixoperator=（constMatrix//destructor~Matrix（）{destroy（）;};//miscellaneousprocedures//returnthenumberofrowintRow（）{returnrow;}//returnthenumberofcolumnintColumn（）{returncolumn;}//returnthereferenceofM[m][n].Thiskindofgrammarislikematlab.Toperator（）（intm,intn）;//overloadtheostreamtemplatefriendstd:

:

ostreamoperator（std:

:

ostreamos,Matrixh）;};

templatevoidMatrix:

:

destroy（）{if（mat!

=NULL）{for（inti=0;i

templateMatrix:

:

Matrix（intn）:

mat（NULL）,row（0）,column（0）{if（n=0）std:

:

cerr\

:

endl;else{mat="new"t*[n];for（inti="0;i

:

endl;>

for（inti=0;i

="j）"mat[i][j]="0;"row="n;"column="n;"}

templateMatrix:

:

Matrix（intm,intn,Tdef）:

mat（NULL）,row（0）,column（0）{if（m0n0）{row=m;column=n;mat=newT*[m];for（inti=0;i

templateMatrix:

:

Matrix（constMatrixM）:

mat（NULL）,row（0）,column（0）{if（M.mat!

=NULL）{mat=newT*[M.row];for（inti=0;i

templateMatrixMatrix:

:

operator=（constMatrixM）{if（this!

=M）{destroy（）;if（M.mat!

=NULL）{mat=newT*[M.row];for（inti=0;i

column=0;}}retur