牛顿差值多项式

实验题目：Newton插值多项式相关知识：

通过n+1个节点的次数不超过n的Newton插值多项式为：

 

x	0	10	20	30	40	50	60	70	80	90	100	110	120
y	5	1	7.5	3	4.5	8.8	15.5	6.5	-5	-10	-2	4.5	7

#include <iostream>
#include <vector>

using namespace std;

//用来保存多项式系数
vector<double> vc; 

// 牛顿差值多项式 数据

double data[][2] = {
	0,5,
	10,1,
	20,7.5,
	30,3,
	40,4.5,
	50,8.8,
    60,15.5,
	70,6.5,
	80,-5,
	90,-10,
	100,-2,
    110,4.5,
	120,7};

//计算出差值函数
void BuildNewtonFunc(double data[][2],const int count)
{
	//计算出每个 基函数 的系数 下面用 小写字母 c+数字 来注释表示
	vc.push_back(data[0][1]); //第一个基函数的系数 为y0 ,同时基函数也是常指函数

	double c0 = vc[0];
	int k ,i;
	double d;
	double u;

	for(k=1;k< count;k++){
		d = data[k][0]-data[k-1][0];  //利用 
		u = vc[k-1];
		for(i=k-2;i>= 0;i--)  //这个for循环利用的是嵌套乘法 或者 霍纳算法 进行分解
		{                     // 目的是用已知的c(i）来推算出 未知的 c(i+i)
			u = u *(data[k][0]-data[i][0])+vc[i];
			d = d *(data[k][0]-data[i][0]);
		}
		vc.push_back((data[k][1]-u)/d);
	}

}
void NewtonFunc(vector<double> vc ,const double x,double &y)
{
	y = 0;
	int i;
	int j;
	double temp;
	for(i=0;i<vc.size();++i){
		temp = vc[i];
		for(j=0;j<i;j++){
			temp *= x - data[j][0];	
		}
		y += temp;
	}
}

int main()
{
	BuildNewtonFunc(data,sizeof(data)/sizeof(data[0]));
	double y;
	NewtonFunc(vc,65,y);
	//cout << y << endl;
	printf("%f",y);


	return 0;
}

 

 

 

 

 

试验要求：利用Newton插值多项式求被插值函数f(x)在点x=65处的近似值。建议：画出Newton插值多项式的曲线。

下面的算法我借鉴的机械工业出版社的一本数值计算的书《数值分析》/.(美) David Kincaid,(美)Ward Cheney著/.王国荣, 俞耀明, 徐兆亮译

 

WIN32代码

 

// application.
//

#include  <windows.h>
#include <iostream>
#include <vector>

LRESULT CALLBACK MainWndProc(HWND hwnd ,UINT message, WPARAM wParam,
							 LPARAM lParam);

//---------------------------------------------------------------------------

using namespace std;

//用来保存多项式系数
vector<double> vc; 

// 牛顿差值多项式 数据

double data[][2] = {
	0,5,
	10,1,
	20,7.5,
	30,3,
	40,4.5,
	50,8.8,
    60,15.5,
	70,6.5,
	80,-5,
	90,-10,
	100,-2,
    110,4.5,
	120,7};

//计算出差值函数
void BuildNewtonFunc(double data[][2],const int count)
{
	//计算出每个 基函数 的系数 下面用 小写字母 c+数字 来注释表示
	vc.push_back(data[0][1]); //第一个基函数的系数 为y0 ,同时基函数也是常指函数

	double c0 = vc[0];
	int k ,i;
	double d;
	double u;

	for(k=1;k< count;k++){
		d = data[k][0]-data[k-1][0];  //利用 
		u = vc[k-1];
		for(i=k-2;i>= 0;i--)  //这个for循环利用的是嵌套乘法 或者 霍纳算法 进行分解
		{                     // 目的是用已知的c(i）来推算出 未知的 c(i+i)
			u = u *(data[k][0]-data[i][0])+vc[i];
			d = d *(data[k][0]-data[i][0]);
		}
		vc.push_back((data[k][1]-u)/d);
	}

}
void NewtonFunc(vector<double> vc ,const double x,double &y)
{
	y = 0;
	int i;
	int j;
	double temp;
	for(i=0;i<vc.size();++i){
		temp = vc[i];
		for(j=0;j<i;j++){
			temp *= x - data[j][0];	
		}
		y += temp;
	}
}


//---------------------------------------------------------------------------





int APIENTRY WinMain(HINSTANCE hInstance,
                     HINSTANCE hPrevInstance,
                     LPSTR     lpCmdLine,
                     int       nCmdShow)
{
	char szClassName[] = "MainWClass";	
	WNDCLASSEX wndclass;
	// 用描述主窗口的参数填充WNDCLASSEX结构
	wndclass.cbSize = sizeof(wndclass); // 结构的大小
	wndclass.style = CS_HREDRAW|CS_VREDRAW;//指定如果大小改变就重画
	wndclass.lpfnWndProc = MainWndProc; // 窗口函数指针
	wndclass.cbClsExtra = 0;//
	wndclass.cbWndExtra = 0; //
	wndclass.hInstance = hInstance;//实例句柄
	wndclass.hIcon = ::LoadIcon(NULL,IDI_APPLICATION);//使用预定义的图标 
	wndclass.hCursor = ::LoadCursor(NULL,IDC_ARROW);//使用预定义的光标 
	wndclass.hbrBackground = (HBRUSH)::GetStockObject(WHITE_BRUSH);
	                                  //使用白色背景画刷
	wndclass.lpszMenuName = NULL;//不指定菜单 
	wndclass.lpszClassName = szClassName;//窗口类的名称 
	wndclass.hIconSm = NULL;//没有类的小图标 
	//注册窗口类
	::RegisterClassEx(&wndclass);
	HWND hwnd = ::CreateWindowEx(
		0,//扩展样式
		szClassName,//类名
		"My First Window",//标题 
		WS_OVERLAPPEDWINDOW,//窗口风格 overlapped window
		CW_USEDEFAULT,
		CW_USEDEFAULT,
		CW_USEDEFAULT,
		CW_USEDEFAULT,
		NULL,//父窗口句柄
		NULL,// 菜单句柄
		hInstance , // 程序实例句柄 
		NULL
		);
	if (hwnd == NULL)
	{	
		::MessageBox(NULL,"ERror","error",MB_OK);
		return -1;
	}
	//显示窗口，刷新客户区
	::ShowWindow(hwnd,nCmdShow);
	::UpdateWindow(hwnd);
	MSG msg ;
	while(::GetMessage(&msg,NULL,0,0))
	{
		::TranslateMessage(&msg);
		::DispatchMessage(&msg);
	}

	return msg.wParam;
}
void DrawLine(HDC hdc,int ax,int ay,int bx,int by)
{
	::MoveToEx(hdc, ax, ay, NULL);
	::LineTo(hdc, bx ,by);
}
LRESULT CALLBACK MainWndProc(HWND hwnd ,UINT message, WPARAM wParam,
							 LPARAM lParam)
{
	char szText[] = "最简单的窗口";
	int cx,cy;
	double i,y,y2;
	switch(message)
	{
	case WM_PAINT:
		{
			HDC hdc;
			PAINTSTRUCT ps;
			hdc = ::BeginPaint(hwnd ,&ps);
			RECT rt;
			::GetClientRect(hwnd,&rt);
	        cx = rt.right;
	        cy = rt.bottom;


			//::TextOut(hdc ,10,10,szText,strlen(szText));

			BuildNewtonFunc(data,sizeof(data)/sizeof(data[0]));
			/*double y;
			NewtonFunc(vc,65,y);*/

			::SetMapMode(hdc,MM_ANISOTROPIC); //设置坐标模式
			::SetWindowExtEx(hdc,600,200,NULL);//设置逻辑单位的大小为1000*1000;
			::SetViewportExtEx(hdc,cx,-cy,NULL);
			::SetViewportOrgEx(hdc,cx/2,cy/2,NULL);//设置坐标原点

			::TextOut(hdc,0,0,"o",1);
		//画X轴
			
			DrawLine(hdc,-cx,0,cy,0);
			DrawLine(hdc,0,cx,0,-cy);

			MoveToEx(hdc,0,0,NULL);
			for(i=0;i<100;i++)
			{
				NewtonFunc(vc,i,y);
				NewtonFunc(vc,i+1,y2);
				
				DrawLine(hdc,i,y,i+1,y2);
			}

			::EndPaint(hwnd ,&ps);
			return 0 ;
		}
	case WM_DESTROY:
		{
			::PostQuitMessage(0);
			return 0 ;
		}
	}

	return ::DefWindowProc(hwnd ,message,wParam ,lParam );
}

效果图：