实验题目: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 );
}
效果图: