数值分析插值上机报告含代码.docx
- 文档编号:9292022
- 上传时间:2023-02-04
- 格式:DOCX
- 页数:8
- 大小:19.28KB
数值分析插值上机报告含代码.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:
:
pair
//linerliststd:
:
list
//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
{
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);//overloadtheostreamtemplate
:
ostreamoperator(std:
:
ostreamos,Matrix
template
:
destroy(){if(mat!
=NULL){for(inti=0;i
template
:
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;"}
template
:
Matrix(intm,intn,Tdef):
mat(NULL),row(0),column(0){if(m0n0){row=m;column=n;mat=newT*[m];for(inti=0;i
template
:
Matrix(constMatrixM):
mat(NULL),row(0),column(0){if(M.mat!
=NULL){mat=newT*[M.row];for(inti=0;i
template
:
operator=(constMatrixM){if(this!
=M){destroy();if(M.mat!
=NULL){mat=newT*[M.row];for(inti=0;i
column=0;}}retur
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 上机 报告 代码
![提示](https://static.bdocx.com/images/bang_tan.gif)