newtoninterpolate.cpp

这是牛顿插值的算法…… 是用于数值计算的 是很不错的 大家可以借鉴一下

#include "StdAfx.h"
#include <assert.h>
#include "NewtonInterpolate.h"

double fac(int n)
{
	double val = 1;
	if(n > 0)
	{
		while(n)
		{
			val *= n--;
		}
	}
	return val;
}


NewtonInterpolate::NewtonInterpolate(void)
{
	m_nNum = -1;
}

NewtonInterpolate::~NewtonInterpolate(void)
{
}

NewtonInterpolate::NewtonInterpolate(int n, double x[], double y[])
{

	assert(n>=0);
	m_nNum = n;
	m_pdNodeX.SetSize(m_nNum+1);
	m_pdNodeY.SetSize(m_nNum+1);
	m_pdDeviationDiagonal.SetSize(m_nNum+1);
	m_pdDeviationLstLine.SetSize(m_nNum+1);

	for(int i =0; i<=n; i++)
	{
		m_pdNodeX[i] = x[i];
		m_pdNodeY[i] = y[i];
	}
	CalculateDeviation();
}

void NewtonInterpolate::CalculateDeviation(void)
{
	int row, col;
	double temp1,temp2;
	m_pdDeviationDiagonal[0] = m_pdDeviationLstLine[0] = m_pdNodeY[0];
	printf("%12.8f\n", m_pdDeviationLstLine[0]);

	for(row=1; row<=m_nNum; row++)
	{
		temp1 = m_pdNodeY[row-1];
		m_pdDeviationLstLine[0] =  m_pdNodeY[row];
//		printf("%12.8f", m_pdNodeY[row]);
		for(col=1; col<= row; col++)
		{
			temp2 = m_pdDeviationLstLine[col];
			m_pdDeviationLstLine[col] = (m_pdDeviationLstLine[col-1]-temp1)/(m_pdNodeX[row]-m_pdNodeX[row-col]);
			temp1 = temp2;
//			printf("%12.8f", m_pdDeviationLstLine[col]);
		}
//		printf("\n");
		m_pdDeviationDiagonal[row] = m_pdDeviationLstLine[row];
	}
}

double NewtonInterpolate::Evaluation(double x)
{
	assert(m_nNum>=0);
	int i = m_nNum;
	double value = m_pdDeviationDiagonal[i];
	while(i-- > 0)
		value = value * (x - m_pdNodeX[i]) + m_pdDeviationDiagonal[i];
	
	return value;	
}

BOOL NewtonInterpolate::AddNode(int degree, double x, double *y)
{
	int row, col;
	double temp1,temp2;
	
	m_nNum += degree;
	int i, j;
	for(i = 0; i<degree; i++)
		m_pdNodeX.Add(x);

	m_pdNodeY.SetSize(m_nNum+1);
	m_pdDeviationDiagonal.SetSize(m_nNum+1);
	m_pdDeviationLstLine.SetSize(m_nNum+1);
	
	row = m_nNum-degree+1;
	for(i = 0; row<=m_nNum; row++, i++)
	{
		if(i == degree-1)
			for(j = 0; j<i; j++)
				m_pdDeviationLstLine[j] = *(y+j)/fac(j);		
		temp1 = m_pdDeviationLstLine[i];
		m_pdDeviationLstLine[i]= *(y+i)/fac(i);	
		
		for(col=i+1; col<= row; col++)
		{
			temp2 = m_pdDeviationLstLine[col];
			m_pdDeviationLstLine[col] = (m_pdDeviationLstLine[col-1]-temp1)/(m_pdNodeX[row]-m_pdNodeX[row-col]);
			temp1 = temp2;
		}
		m_pdDeviationDiagonal[row] = m_pdDeviationLstLine[row];
	}

	return TRUE;//参数不正确时返回FALSE, 待开发
}

NewtonInterpolate::NewtonInterpolate(int n, int degree[], double x[], double *y[])
{
	m_nNum = -1;
	int i;
	for(i = 0; i <= n; i++)
		AddNode(degree[i], x[i], y[i]);
}

void NewtonInterpolate::PrintPolynomial()
{
	printf("Newton插值多项式的系数为:\n");
	for(int i = 0; i <= m_nNum; i++)
		printf("%8.4f", m_pdDeviationDiagonal[i]);

	printf("\n");
}