interpolationdoc.cpp

刚学数值分析

// InterpolationDoc.cpp : implementation of the CInterpolationDoc class
//

#include "stdafx.h"
#include "Interpolation.h"

#include "InterpolationDoc.h"

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

/////////////////////////////////////////////////////////////////////////////
// CInterpolationDoc

IMPLEMENT_DYNCREATE(CInterpolationDoc, CDocument)

BEGIN_MESSAGE_MAP(CInterpolationDoc, CDocument)
	//{{AFX_MSG_MAP(CInterpolationDoc)
		// 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
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CInterpolationDoc construction/destruction

CInterpolationDoc::CInterpolationDoc()
{
	// TODO: add one-time construction code here
	m_nNumHits	=	10;
	m_nNumDraw	=	1000;

	m_dblXmin	=	-1.0;
	m_dblXmax	=	1.0;
	m_dblh		=	(m_dblXmax-m_dblXmin)/m_nNumHits;

	m_dblX	=	new double[m_nNumDraw+1];
	double	dblStep	=	(m_dblXmax-m_dblXmin)/m_nNumDraw;
	for(int i=0; i<=m_nNumDraw; i++)
	{
		m_dblX[i] = m_dblXmin + i*dblStep;
	}

	m_dblF	=	new double[m_nNumDraw+1];
	m_dblL	=	new double[m_nNumDraw+1];
	m_dblS	=	new double[m_nNumDraw+1];

	m_x	=	new double[m_nNumHits+1];
	m_y	=	new double[m_nNumHits+1];
	for(i=0; i<=m_nNumHits; i++)
	{
		m_x[i] = m_dblXmin + i*m_dblh;
		m_y[i] = FuncF(m_x[i]);
	}
}

CInterpolationDoc::~CInterpolationDoc()
{
}

BOOL CInterpolationDoc::OnNewDocument()
{
	if (!CDocument::OnNewDocument())
		return FALSE;

	// TODO: add reinitialization code here
	// (SDI documents will reuse this document)
	GetF();
	GetL();
	GetS();

	return TRUE;
}



/////////////////////////////////////////////////////////////////////////////
// CInterpolationDoc serialization

void CInterpolationDoc::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		// TODO: add storing code here
	}
	else
	{
		// TODO: add loading code here
	}
}

/////////////////////////////////////////////////////////////////////////////
// CInterpolationDoc diagnostics

#ifdef _DEBUG
void CInterpolationDoc::AssertValid() const
{
	CDocument::AssertValid();
}

void CInterpolationDoc::Dump(CDumpContext& dc) const
{
	CDocument::Dump(dc);
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CInterpolationDoc commands
void CInterpolationDoc::GetF()
{
	for(int i=0; i<=m_nNumDraw; i++)
	{
		m_dblF[i] = FuncF(m_dblX[i]);
	}
}

void CInterpolationDoc::GetL()
{
	for(int i=0; i<=m_nNumDraw; i++)
	{
		m_dblL[i] = FuncL(m_nNumHits, m_dblX[i]);
	}
}

void CInterpolationDoc::GetS()
{
	double*	Mu		=	new double[m_nNumHits+1];
	double*	LambDa	=	new double[m_nNumHits+1];
	double*	d		=	new double[m_nNumHits+1];
	double*	B		=	new double[m_nNumHits+1];


	for(int i=0; i<=m_nNumHits; i++)
	{
		Mu[i]	=	0.5;
	}
	Mu[m_nNumHits] = 0.0;

	for(i=0; i<=m_nNumHits; i++)
	{
		LambDa[i]	=	0.5;
	}
	LambDa[0] = 0.0;

	d[0] = d[m_nNumHits] = 0.0;
	for(i=1; i<=m_nNumHits-1; i++)
	{
		d[i] = 3*(m_y[i+1] + m_y[i-1] - 2*m_y[i])/(m_dblh*m_dblh);
	}

	for(i=0; i<=m_nNumHits; i++)
	{
		B[i] = 2.0;
	}

	ThreeDiag TD(m_nNumHits, Mu, B, LambDa, d);
	double* M = TD.GetSolution();


	for(i=0; i<=m_nNumDraw; i++)
	{
		for(int j=0; j<=m_nNumHits-1; j++)
		{
			if(m_x[j]<=m_dblX[i] && m_dblX[i]<=m_x[j+1])
			{
				double xxx = 0;
				xxx += M[j]*(m_x[j+1]-m_dblX[i])*(m_x[j+1]-m_dblX[i])*(m_x[j+1]-m_dblX[i]);
				xxx += M[j+1]*(m_dblX[i]-m_x[j])*(m_dblX[i]-m_x[j])*(m_dblX[i]-m_x[j]);
				xxx += (6*m_y[j]-M[j]*m_dblh*m_dblh)*(m_x[j+1]-m_dblX[i]);
				xxx += (6*m_y[j+1]-M[j+1]*m_dblh*m_dblh)*(m_dblX[i]-m_x[j]);
				m_dblS[i] = xxx/(6*m_dblh);
			}
		}
	}
}

double CInterpolationDoc::FuncF(double x)
{
	return 1.0/(1.0+9.0*x*x);
}

double CInterpolationDoc::FuncL(int n, double x)
{
	double sum = 0;
	for(int i=0; i<=n; i++)
	{
		double sum2 = 1;
		for(int j=0; j<=m_nNumHits; j++)
		{
			if(j != i)
			{
				double xxx = (x-m_x[j])/(m_x[i]-m_x[j]);
				sum2 *= xxx;
			}
		}

		sum += m_y[i]*sum2;
	}

	return sum;
}

double CInterpolationDoc::FuncS(int n, double x, double* M)
{
	for(int j=0; j<=m_nNumHits-1; j++)
	{
		if(m_x[j]<=x && x<=m_x[j+1])
		{
			double xxx = 0;
			xxx += M[j]*(m_x[j+1]-x)*(m_x[j+1]-x)*(m_x[j+1]-x);
			xxx += M[j+1]*(x-m_x[j])*(x-m_x[j])*(x-m_x[j]);
			xxx += (6*m_y[j]-M[j]*m_dblh*m_dblh)*(m_x[j+1]-x);
			xxx += (6*m_y[j+1]-M[j+1]*m_dblh*m_dblh)*(x-m_x[j]);
			xxx = xxx/(6*m_dblh);

//			xxx += M[j]*(m_x[j+1]-x)*(m_x[j+1]-x)*(m_x[j+1]-x)/(6*m_dblh);
//			xxx += M[j+1]*(x-m_x[j])*(x-m_x[j])*(x-m_x[j])/(6*m_dblh);
//			xxx += (m_y[j]-M[j]*m_dblh*m_dblh/6)*(m_x[j+1]-x)/m_dblh;
//			xxx += (m_y[j+1]-M[j+1]*m_dblh*m_dblh/6)*(x-m_x[j])/m_dblh;
			return	 xxx;
		}
	}

	return 0;
}