threesplineview.cpp

数值分析中三次样条插值算法

// ThreeSplineView.cpp : implementation of the CThreeSplineView class
//

#include "stdafx.h"
#include "ThreeSpline.h"
#include "math.h"
#include "ThreeSplineDoc.h"
#include "ThreeSplineView.h"
#define N 6
#define  M 5*(N-1)+1
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CThreeSplineView

IMPLEMENT_DYNCREATE(CThreeSplineView, CView)

BEGIN_MESSAGE_MAP(CThreeSplineView, CView)
	//{{AFX_MSG_MAP(CThreeSplineView)
		// NOTE - the ClassWizard will add and remove mapping macros here.
		//    DO NOT EDIT what you see in these blocks of generated code!
	//}}AFX_MSG_MAP
	// Standard printing commands
	ON_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()

/////////////////////////////////////////////////////////////////////////////
// CThreeSplineView construction/destruction

CThreeSplineView::CThreeSplineView()
{
	// TODO: add construction code here

}

CThreeSplineView::~CThreeSplineView()
{
}

BOOL CThreeSplineView::PreCreateWindow(CREATESTRUCT& cs)
{
	// TODO: Modify the Window class or styles here by modifying
	//  the CREATESTRUCT cs

	return CView::PreCreateWindow(cs);
}

/////////////////////////////////////////////////////////////////////////////
// CThreeSplineView drawing

void CThreeSplineView::OnDraw(CDC* pDC)
{
	CThreeSplineDoc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	// TODO: add draw code for native data here
	//三次样条插值函数，选取函数y=pow(x,4)-10*pow(x,3)，
	//满足在最前和最后两端点处二次导数为0，其中scaleX[N]用于放插值节点的X值，
	//y[N]用于存放插值节点对应的X值
	//initX[M],initY[M],initZ[M]分别放画图时的横坐标和纵坐标.initY[M]放插值后的函数值,
	//initZ[M]放插值前的函数值.
	float initX[M];
    float initY[M];
    float initZ[M];
    float  scaleX[N];
	float x[N];
    float y[N];
	float alpha[N];
	float beta[N];
	float m[N];
	float A[N];
    float B[N];
	float h[N];
	int i,j,k;
	
    for(i=0;i<N;i++)
	{
		x[i]=i;
		scaleX[i]=i;
		y[i]=pow(x[i],4)-10*pow(x[i],3);
	}
	
	for (i=0;i<N-1;i++)
	{
		h[i]=scaleX[i+1]-scaleX[i];
	}
    alpha[0]=1;
	alpha[N-1]=0;
    beta[0]=3*(y[1]-y[0])/h[0];
    beta[N-1]=3*(y[N-1]-y[N-2])/h[N-2];
	for(i=1;i<N-1;i++)
	{
		alpha[i]=h[i-1]/(h[i-1]+h[i]);
		beta[i]=3*((1-alpha[i])*(y[i]-y[i-1])/h[i-1]+alpha[i]*(y[i+1]-y[i])/h[i]);
	}
	A[0]=-alpha[0]/2;
	B[0]=beta[0]/2;
	for(i=1;i<N;i++)
	{
       A[i]=-alpha[i]/(2+(1-alpha[i])*A[i-1]);
	   B[i]=(beta[i]-(1-alpha[i])*B[i-1])/(2+(1-alpha[i])*A[i-1]);
	}
	m[N-1]=(beta[N-1]-(1-alpha[N-1])*B[N-2])/(2+(1-alpha[N-1])*A[N-2]);
	for(i=N-2;i>=0;i--)
	{
		m[i]=A[i]*m[i+1]+B[i];
	}
	for (i=0;i<M;i++)
	{
		initX[i]=float(i)*0.2;
	}
	for(i=0;i<M;i++)
	{
		j=0;
		while((initX[i]>=scaleX[j]&&initX[i]>scaleX[j+1]))
		{
			j++;
		}
		initY[i]=(1+2*(initX[i]-scaleX[j])/h[j])*((initX[i]-scaleX[j+1])/(-h[j]))*((initX[i]-scaleX[j+1])/(-h[j]))*y[j]+
			(1+2*(initX[i]-scaleX[j+1])/(-h[j]))*((initX[i]-scaleX[j])/(h[j]))*((initX[i]-scaleX[j])/(h[j]))*y[j+1]+
            (initX[i]-scaleX[j])*((initX[i]-scaleX[j+1])/(-h[j]))*((initX[i]-scaleX[j+1])/(-h[j]))*m[j]+
    		(initX[i]-scaleX[j+1])*((initX[i]-scaleX[j])/(h[j]))*((initX[i]-scaleX[j])/(h[j]))*m[j+1];
        initZ[i]=pow(initX[i],4)-10*pow(initX[i],3);
	}

	 //分别取得X和Y最大值,最小值

       int  maxX=initX[0];

       int  minX=initX[0];

       float maxY=initY[0];

       float minY=initY[0];
	   float maxZ=initZ[0];

       float minZ=initZ[0];

       for( i=0;i<M;i++)

       {

              if(initX[i]>maxX)

                     maxX=initX[i]; 

              if(initX[i]<minX)

                     minX=initX[i];

              if(initY[i]>maxY)

                     maxY=initY[i];

              if(initY[i]<minY)

                     minY=initY[i];
			   if(initZ[i]>maxZ)

                     maxZ=initZ[i];

              if(initZ[i]<minZ)

                     minZ=initZ[i];
             
       }

      //如果原点必须在X轴上，加上下面2行，否则注释掉

       if(minX>0)

              minX=0;

      //如果原点必须在Y轴上，加上下面2行，否则注释掉

       if(minY>0)

              minY=0;
         
       if(minZ>0)

              minZ=0;
	  
 

       //确定图象显示大小

       int width=800;
       int height=300;
	   //要显示两个图象,只要把下面一个图象其高度都增加一个值即可.

       //确定坐标图四周预留的空白大小
       const int mytop=30;
	   const int mybottom=40;
	   const int myleft=80;
       const int myright=50;

       //确定X，Y,Z轴每单位显示宽度
       float intervalX=(width-myleft-myright)/(maxX-minX);
       float intervalY=(height-mybottom-mytop)/(maxY-minY);
       float intervalZ=(height-mybottom-mytop)/(maxZ-minZ);
        
       //绘制曲线。由于绘图坐标的Y轴是向下延升，所以这里每个点的Y值是用
       //图象高度减去y值大小。
       pDC->MoveTo(int(myleft+(initX[0]-minX)*intervalX), 
             int(height-(mybottom+(initY[0]-minY)*intervalY)));
       for(i=0;i<M;i++)
	   {

       pDC->LineTo(int(myleft+(initX[i]-minX)*intervalX), 
		   int(height-(mybottom+(initY[i]-minY)*intervalY)));//画线,与moveto连用    
       }

       //绘制X轴,从图形区域最左端到最右端
       float bottomY=0;
	   float leftX=0;
       float bottomZ=0;

       //bottomY表示X轴的y值，leftX表示Y轴的x值
       if(minY>0)
           bottomY=minY;
       if(minX>0)
           leftX=minX;
	   if(minZ>0)
           bottomZ=minZ;
       pDC->MoveTo(int(myleft),int(height-(mybottom+(bottomY-minY)*intervalY)));
       pDC->LineTo(int(width-myright),int(height-(mybottom+(bottomY-minY)*intervalY)));

       //Y轴从图形区域最底端到最顶端
       pDC->MoveTo(int(myleft+(leftX-minX)*intervalX),int(height-mybottom));
       pDC->LineTo(int(myleft+(leftX-minX)*intervalX),int(mytop));

       //确定显示刻度个数
       const int count=5;
	   int countx=M-1;

       //确定每个显示刻度之间的宽度
       float spaceX=(width-myleft-myright)/countx;
       float spaceY=(height-mybottom-mytop)/count;
       float spaceZ=(height-mybottom-mytop)/count;

       //绘制刻度和刻度值
       CString str;

       //X轴
      for(i=0;i<=countx;i++)//yhy
	  {
              //str.Format("%d",minX+i*(maxX-minX)/countx);
              str.Format("%3.1f",i*0.2);
              pDC->MoveTo(int(myleft+spaceX*i),int(height-(mybottom+(bottomY-minY)*intervalY)));

              pDC->LineTo(int(myleft+spaceX*i),int(height-(mybottom+(bottomY-minY)*intervalY+5)));

              pDC->TextOut(int(myleft+spaceX*i-10),

                     int(height-(mybottom+(bottomY-minY)*intervalY-5)),str);
       }
        //Y轴

        for(i=0;i<=count;i++)

       {

             str.Format("%.1f",minY+i*(maxY-minY)/count);

              pDC->MoveTo(int(myleft+(leftX-minX)*intervalX),int(height-(mybottom+spaceY*i)));

              pDC->LineTo(int(myleft+(leftX-minX)*intervalX+5),int(height-(mybottom+spaceY*i)));

              pDC->TextOut(int(myleft+(leftX-minX)*intervalX-30),
//
                  int(height-(mybottom+spaceY*i+8)),str);

       }

      //绘制X，Y轴的变量名

       pDC->TextOut(width/3,height-20,"图1：三次样条插值函数");

       pDC->TextOut(0,height/2,"幅值");
		 //X轴
//test start
		pDC->MoveTo(int(myleft+(initX[0]-minX)*intervalX), 

                     int(height+300-(mybottom+(initZ[0]-minZ)*intervalZ)));

       for(i=0;i<M;i++)

       {

            pDC->LineTo(int(myleft+(initX[i]-minX)*intervalX), 

                    int(height+300-(mybottom+(initZ[i]-minZ)*intervalZ)));
       }
      //绘制X轴和Y轴
       pDC->MoveTo(int(myleft),int(height+300-(mybottom+(bottomZ-minZ)*intervalZ)));

       pDC->LineTo(int(width-myright),int(height+300-(mybottom+(bottomZ-minZ)*intervalZ)));

       //Y轴从图形区域最底端到最顶端

       pDC->MoveTo(int(myleft+(leftX-minX)*intervalX),int(height+300-mybottom));

       pDC->LineTo(int(myleft+(leftX-minX)*intervalX),int(mytop+300));

       for(i=0;i<=countx;i++)

       {

              //str.Format("%d",minX+i*(maxX-minX)/countx);
		      str.Format("%3.1f",i*0.2);

              pDC->MoveTo(int(myleft+spaceX*i),int(height+300-(mybottom+(bottomZ-minZ)*intervalZ)));

              pDC->LineTo(int(myleft+spaceX*i),int(height+300-(mybottom+(bottomZ-minZ)*intervalZ+5)));//每个MOVETO ,LINeto 画刻度上的小竖线

              pDC->TextOut(int(myleft+spaceX*i-10),

                     int(height+300-(mybottom+(bottomZ-minZ)*intervalZ-5)),str);

       }
      
 

       //Y轴

        for(i=0;i<=count;i++)

       {

              str.Format("%.1f",minZ+i*(maxZ-minZ)/count);

              pDC->MoveTo(int(myleft+(leftX-minX)*intervalX),int(height+300-(mybottom+spaceZ*i)));

              pDC->LineTo(int(myleft+(leftX-minX)*intervalX+5),int(height+300-(mybottom+spaceZ*i)));

              pDC->TextOut(int(myleft+(leftX-minX)*intervalZ-30),
//
                  int(height+300-(mybottom+spaceZ*i+8)),str);

       }

          //绘制X，Y轴的变量名

       pDC->TextOut(width/3,height+300-20,"图2：原函数(y=pow(x,4)-10*pow(x,3))");

       pDC->TextOut(0,height/2+300,"幅值");
}

/////////////////////////////////////////////////////////////////////////////
// CThreeSplineView printing

BOOL CThreeSplineView::OnPreparePrinting(CPrintInfo* pInfo)
{
	// default preparation
	return DoPreparePrinting(pInfo);
}

void CThreeSplineView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: add extra initialization before printing
}

void CThreeSplineView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: add cleanup after printing
}

/////////////////////////////////////////////////////////////////////////////
// CThreeSplineView diagnostics

#ifdef _DEBUG
void CThreeSplineView::AssertValid() const
{
	CView::AssertValid();
}

void CThreeSplineView::Dump(CDumpContext& dc) const
{
	CView::Dump(dc);
}

CThreeSplineDoc* CThreeSplineView::GetDocument() // non-debug version is inline
{
	ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CThreeSplineDoc)));
	return (CThreeSplineDoc*)m_pDocument;
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CThreeSplineView message handlers