当前位置:文档之家› 实验六:分段三次Hermit插值

实验六:分段三次Hermit插值

《数值分析》实验报告实验序号:实验六题目名称: 分段三次Hermit插值学号: 姓名: 、任课教师: 马季骕专业班级:计算机科学与技术(非师范)1、题目描述这是一种光滑的分段插值,给定[a,b]上的节点a=x0<x1<x2<…xn-1<xn=b和微商值f’(xi)=y’i(i=0,1,2…,n),作一个分段三次Hermit插值函数H(x),要求满足条件:(1)H(xi)=yi,H’(xi)=yi(i=0,1,2,…,n)。

(2)在每一个小区间[xi,xi+1]上是三次多项式。

根据前的结果做出个点的插值基函数2、算法分析●分段三次Hermit插值的算法思想:分段三次Hermit插值的做法是在每一个小区间上作三次Hermit插值,因此在每一个插值节点上都需要构造两个插值基函数,然后再作它们的线性组合。

分段三次Hermit插值基函数如下:H(x)=Σ(yihi(x)+y’iHi(x))●具体程序设计:for(i=0;i<=10;i++){ax[i]=-1+((2/10.0)*i); //选取节点ay[i]=1.0/(1+25*ax[i]*ax[i]); //函数值a[i]=(-50*ax[i])/((1+25*ax[i]*ax[i])*(1+25*ax[i]*ax[i]));//一阶导数值}x1=-1;y1=1/(1+25*x1*x1);while(x1<=1){m=0;x1=x1+0.00001;for(i=0;i<=10;i++){if(i==0&&x1>=ax[0]&&x1<=ax[1]){n=(1+2*((x1-ax[0])/(ax[1]-ax[0])))*((x1-ax[1])/(ax[0]-ax[ 1]))*((x1-ax[1])/(ax[0]-ax[1]));// 插值基函数h(x)tt=ay[0]*n; //插值基函数与其函数值线性组合h=(x1-ax[0])*((x1-ax[1])/(ax[0]-ax[1]))*((x1-ax[1])/(ax[0 ]-ax[1]));//插值基函数H(x)tt1=a[0]*h;//插值基函数与其一阶导数值线性组合m=tt+tt1+m;//把插值基函数做线性组合}else{if(x1>=ax[i-1]&&x1<=ax[i]){n=(1+2*((x1-ax[i])/(ax[i-1]-ax[i])))*((x1-ax[i-1])/(ax[i]-ax[i-1]))*((x1-ax[i-1])/(ax[i]-ax[i-1]));tt=ay[i]*n;h=(x1-ax[i])*((x1-ax[i-1])/(ax[i]-ax[i-1]))*((x1-ax[i -1])/(ax[i]-ax[i-1]));tt1=a[i]*h;m=tt+tt1+m;}else{if(x1>ax[i]&&x1<=ax[i+1]){n=(1+2*((x1-ax[i])/(ax[i+1]-ax[i])))*((x1-ax[i+ 1])/(ax[i]-ax[i+1]))*((x1-ax[i+1])/(ax[i]-ax[i+1]));tt=ay[i]*n;h=(x1-ax[i])*((x1-ax[i+1])/(ax[i]-ax[i+1]))*((x1 -ax[i+1])/(ax[i]-ax[i+1]));tt1=a[i]*h;m=tt+tt1+m;}else{if(i==9&&x1>ax[9]&&x1<=ax[10]){n=(1+2*((x1-ax[10])/(ax[9]-ax[10])))*((x1-ax [9])/(ax[10]-ax[9]))*((x1-ax[9])/(ax[10]-ax[9]));tt=ay[i]*n;h=(x1-ax[10])*((x1-ax[9])/(ax[10]-ax[9]))*(( x1-ax[9])/(ax[10]-ax[9]));tt1=a[i]*h;m=tt+tt1+m;}}}}}}3、运行结果截图:4、程序代码// SHIYAN456View.cpp : implementation of the CSHIYAN456View class//#include "stdafx.h"#include "SHIYAN456.h"#include "SHIYAN456Doc.h"#include "SHIYAN456View.h"#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif///////////////////////////////////////////////////////////////// ////////////// CSHIYAN456ViewIMPLEMENT_DYNCREATE(CSHIYAN456View, CView)BEGIN_MESSAGE_MAP(CSHIYAN456View, CView)//{{AFX_MSG_MAP(CSHIYAN456View)ON_COMMAND(ID_FFunction, OnFFunction)ON_COMMAND(ID_Lagrange, OnLagrange)ON_COMMAND(ID_Subsection, OnSubsection)ON_COMMAND(ID_Hermite, OnHermite)//}}AFX_MSG_MAP// Standard printing commandsON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview) END_MESSAGE_MAP()///////////////////////////////////////////////////////////////// ////////////// CSHIYAN456View construction/destructionCSHIYAN456View::CSHIYAN456View(){// TODO: add construction code here}CSHIYAN456View::~CSHIYAN456View(){}BOOL CSHIYAN456View::PreCreateWindow(CREATESTRUCT& cs){// TODO: Modify the Window class or styles here by modifying // the CREATESTRUCT csreturn CView::PreCreateWindow(cs);}///////////////////////////////////////////////////////////////// ////////////// CSHIYAN456View drawingvoid CSHIYAN456View::OnDraw(CDC* pDC){CSHIYAN456Doc* pDoc = GetDocument();ASSERT_VALID(pDoc);// TODO: add draw code for native data hereAfxGetMainWnd()->SetWindowText("实验四五六函数图像");for(int k=650;k<=690;k++){pDC->SetPixel(k,55,RGB(255,0,0));pDC->SetPixel(k,85,RGB(0,255,0));pDC->SetPixel(k,115,RGB(0,0,255));pDC->SetPixel(k,145,RGB(0,255,255));}pDC->TextOut(700,50,"原函数图像");pDC->TextOut(700,80,"Lagrange插值函数图像");pDC->TextOut(700,110,"分段线性插值函数图像");pDC->TextOut(700,140,"Hermite插值函数图像");for(int i=6;i<=600;i++)pDC->SetPixel(400,i,RGB(0,0,0));pDC->TextOut(395,4,"y");for(int j=100;j<=700;j++)pDC->SetPixel(j,500,RGB(0,0,0));pDC->TextOut(700,500,"x");//pDC->MoveTo(400,500);}///////////////////////////////////////////////////////////////// ////////////// CSHIYAN456View printingBOOL CSHIYAN456View::OnPreparePrinting(CPrintInfo* pInfo){// default preparationreturn DoPreparePrinting(pInfo);}void CSHIYAN456View::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/){// TODO: add extra initialization before printing}void CSHIYAN456View::OnEndPrinting(CDC* /*pDC*/, CPrintInfo*/*pInfo*/){// TODO: add cleanup after printing}///////////////////////////////////////////////////////////////// ////////////// CSHIYAN456View diagnostics#ifdef _DEBUGvoid CSHIYAN456View::AssertValid() const{CView::AssertValid();}void CSHIYAN456View::Dump(CDumpContext& dc) const{CView::Dump(dc);}CSHIYAN456Doc* CSHIYAN456View::GetDocument() // non-debug version is inline{ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CSHIYAN456Doc)));return (CSHIYAN456Doc*)m_pDocument;}#endif //_DEBUG///////////////////////////////////////////////////////////////// ////////////// CSHIYAN456View message handlersvoid CSHIYAN456View::OnFFunction(){// TODO: Add your command handler code hereCClientDC dr(this);COLORREF rgb=RGB(255,0,0);double x1,y1,x,y;x1=-1.0;y1=1/(1+25*x1*x1);x=x1*200+400;y=-y1*200+500;while(x1<=1){dr.MoveTo(int(x),int(y));x1=x1+0.00001;y1=1/(1+25*x1*x1);x=x1*200+400;y=-y1*200+500;dr.SetPixel(int(x),int(y),rgb);}}void CSHIYAN456View::OnLagrange(){// TODO: Add your command handler code hereCClientDC dr(this);COLORREF rgb=RGB(0,255,0);int i,j;double x1=0,y1=0,x=0,y=0,m=0,n=0,ax[100],ay[100]; for(i=0;i<=10;i++){ax[i]=-1+((2/10.0)*i);ay[i]=1.0/(1+25*ax[i]*ax[i]);}x1=-1;y1=1/(1+25*x1*x1);x=x1*200+400;y=-y1*200+500;dr.MoveTo(int(x),int(y));while(x1<=1){m=0;for(i=0;i<=10;i++){n=1;for(j=0;j<=10;j++){if(i!=j)n=((x1-ax[j])/(ax[i]-ax[j]))*n;}m=ay[i]*n+m;}x=x1*200+400;y=-m*200+500;dr.SetPixel(int(x),int(y),rgb);x1=x1+0.00001;}}void CSHIYAN456View::OnSubsection(){// TODO: Add your command handler code hereCClientDC dr(this);COLORREF rgb=RGB(0,0,255);int i;double x1=0,y1=0,x=0,y=0,m=0,n=0,ax[100],ay[100]; for(i=0;i<=20;i++){ax[i]=-1+((2/20.0)*i);ay[i]=1.0/(1+25*ax[i]*ax[i]);}x1=-1;y1=1/(1+25*x1*x1);x=x1*200+400;y=-y1*200+500;dr.MoveTo(int(x),int(y));while(x1<=1){m=0;x1=x1+0.00001;for(i=0;i<=20;i++){if(i==0&&x1>=-1&&x1<=-0.9){n=ay[0]*((x1-ax[1])/(ax[0]-ax[1]));m=n+m;}else{if(x1>=ax[i-1]&&x1<=ax[i]){n=ay[i]*((x1-ax[i-1])/(ax[i]-ax[i-1])); m=n+m;}else{if(x1>ax[i]&&x1<=ax[i+1]){n=ay[i]*((x1-ax[i+1])/(ax[i]-ax[i+1]));m=m+n;}else{if(i==19&&x1>ax[19]&&x1<=ax[20]){n=ay[20]*((x1-ax[19])-(ax[20]-ax[19]) );m=n+m;}}}}}x=x1*200+400;y=-m*200+500;dr.SetPixel(int(x),int(y),rgb);}}void CSHIYAN456View::OnHermite(){// TODO: Add your command handler code hereCClientDC dr(this);COLORREF rgb=RGB(0,255,255);int i;doublex1=0,y1=0,x=0,y=0,m=0,n=0,tt=0,tt1=0,h=0,ax[100],ay[100],a[100];for(i=0;i<=10;i++){ax[i]=-1+((2/10.0)*i);ay[i]=1.0/(1+25*ax[i]*ax[i]);a[i]=(-50*ax[i])/((1+25*ax[i]*ax[i])*(1+25*ax[i]*ax[i]));}x1=-1;y1=1/(1+25*x1*x1);x=x1*200+400;y=-y1*200+500;dr.MoveTo(int(x),int(y));while(x1<=1){m=0;x1=x1+0.00001;for(i=0;i<=10;i++){if(i==0&&x1>=ax[0]&&x1<=ax[1]){n=(1+2*((x1-ax[0])/(ax[1]-ax[0])))*((x1-ax[1])/( ax[0]-ax[1]))*((x1-ax[1])/(ax[0]-ax[1]));tt=ay[0]*n;h=(x1-ax[0])*((x1-ax[1])/(ax[0]-ax[1]))*((x1-ax[ 1])/(ax[0]-ax[1]));tt1=a[0]*h;m=tt+tt1+m;}else{if(x1>=ax[i-1]&&x1<=ax[i]){n=(1+2*((x1-ax[i])/(ax[i-1]-ax[i])))*((x1-ax [i-1])/(ax[i]-ax[i-1]))*((x1-ax[i-1])/(ax[i]-ax[i-1]));tt=ay[i]*n;h=(x1-ax[i])*((x1-ax[i-1])/(ax[i]-ax[i-1]))* ((x1-ax[i-1])/(ax[i]-ax[i-1]));tt1=a[i]*h;m=tt+tt1+m;}else{if(x1>ax[i]&&x1<=ax[i+1]){n=(1+2*((x1-ax[i])/(ax[i+1]-ax[i])))*(( x1-ax[i+1])/(ax[i]-ax[i+1]))*((x1-ax[i+1])/(ax[i] -ax[i+1]));tt=ay[i]*n;h=(x1-ax[i])*((x1-ax[i+1])/(ax[i]-ax[i+ 1]))*((x1-ax[i+1])/(ax[i]-ax[i+1]));tt1=a[i]*h;m=tt+tt1+m;}else{if(i==9&&x1>ax[9]&&x1<=ax[10]){n=(1+2*((x1-ax[10])/(ax[9]-ax[10])) )*((x1-ax[9])/(ax[10]-ax[9]))*((x1-ax[9])/(ax[10 ]-ax[9]));tt=ay[i]*n;h=(x1-ax[10])*((x1-ax[9])/(ax[10]-a x[9]))*((x1-ax[9])/(ax[10]-ax[9]));tt1=a[i]*h;m=tt+tt1+m;}}}}}x=x1*200+400;y=-m*200+500;dr.SetPixel(int(x),int(y),rgb);}}5、总结体会:通过本次实验我对分段三次Hermit插值有了更深刻更全面的掌握,它在给定了节点处的函数值和导数值以后,构造了一个整体上具有一阶连续微商的插值函数。

相关主题