牛顿差值多项式

实验题目：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 windowCW_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 );
}

效果图：