lagrange.cpp

计算方法程序常微分方程的数值解法课堂讲义

//////////////////////////////////////////////
//	§2  插值方法
//
//	程序2.1 — Lagrange 插值

////////////////////////////////////////////////////////////////////////
//
//本程序最后将显示函数表和插值函数图形曲线，是TC版本，必须用TC,BC编译.
//如果对Visual C/C++的Windows程序有兴趣，请发e-mail索取,其中包括了所有
//其中包括了所有插值和曲线拟合等功能强大的程序

#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <string.h>
#include <graphics.h>
#include <math.h>

//点结构定义
typedef struct {
	float x;
	float y;
} POINT;

//插值函数原型
float Lagrange(POINT Table[],int n,float x);

//图形初始化
void InitGraph()
{
	int gdriver = DETECT, gmode;

	initgraph(&gdriver, &gmode, "");
}

//画线
void line(int x1,int y1,int x2,int y2,int color)
{
	setcolor(color);
	line(x1,y1,x2,y2);
}

//字符串显示
void outtextxy(int x,int y,char text[],int color)
{
	setcolor(color);
	outtextxy(x,y,text);
}

//实数显示
void outtextxy(int x,int y,float RealNumber,int Length,int color)
{
	char String[20];

	setcolor(color);
	gcvt(RealNumber,Length,String);//实数转换为字符串
	outtextxy(x,y,String);
}

//表格线显示
void ShowTable(POINT Table[],int n)
{
	int i;
	char String_x[20],String_y[20];

	line(0,20,639,20,7);
	line(30,0,30,40);
	outtextxy(10,8,"x",14);
	outtextxy(10,28,"y",14);

	for(i=0;i<n;++i)
	{
		outtextxy(36+i*80,8,Table[i].x,6,14);
		outtextxy(36+i*80,28,Table[i].y,6,14);
	}
}

//得到函数表中x值的绝对值最大值,用来计算坐标的刻度单位
float GetMaximum_x(POINT Table[],int n)
{
	int i;
	int Max_x;

	Max_x=0.0;
	for(i=0;i<n;++i)
		if(Max_x<fabs(Table[i].x)) Max_x=fabs(Table[i].x);
	return(Max_x);
}

//得到函数表中y值的绝对值最大值,用来计算坐标的刻度单位
float GetMaximum_y(POINT Table[],int n)
{
	int i;
	int Max_y;

	Max_y=0.0;
	for(i=0;i<n;++i)
		if(Max_y<fabs(Table[i].y)) Max_y=fabs(Table[i].y);
	return(Max_y);
}

//画直角坐标
void DrawCoordinate(int X1,int Y1,int X2,int Y2,
					float x1,float y1,float x2,float y2)
{
	char String[20];

	line(X1,(Y1+Y2)/2,X2,(Y1+Y2)/2,7);
	line(X2,(Y1+Y2)/2,X2-10,(Y1+Y2)/2-2);
	line(X2,(Y1+Y2)/2,X2-10,(Y1+Y2)/2+2);
	line((X1+X2)/2,Y1,(X1+X2)/2,Y2);
	line((X1+X2)/2,Y1,(X1+X2)/2-2,Y1+10);
	line((X1+X2)/2,Y1,(X1+X2)/2+2,Y1+10);

	line((X1+X2)/2+(X2-X1)/4,(Y1+Y2)/2-4,
		 (X1+X2)/2+(X2-X1)/4,(Y1+Y2)/2);
	gcvt((x1+x2)/2+(x2-x1)/4,6,String);
	outtextxy((X1+X2)/2+(X2-X1)/4-3,(Y1+Y2)/2+4,String);

	line((X1+X2)/2-(X2-X1)/4,(Y1+Y2)/2-4,
		 (X1+X2)/2-(X2-X1)/4,(Y1+Y2)/2);
	gcvt(-((x1+x2)/2+(x2-x1)/4),6,String);
	outtextxy((X1+X2)/2-(X2-X1)/4-3,(Y1+Y2)/2+4,String);

	line((X1+X2)/2+4,(Y1+Y2)/2-(Y2-Y1)/4,
		 (X1+X2)/2,(Y1+Y2)/2-(Y2-Y1)/4);
	gcvt(-((y1+y2)/2+(y2-y1)/4),6,String);
	outtextxy((X1+X2)/2+8,(Y1+Y2)/2-(Y2-Y1)/4-2,String);

	line((X1+X2)/2+4,(Y1+Y2)/2+(Y2-Y1)/4,
		 (X1+X2)/2,(Y1+Y2)/2+(Y2-Y1)/4);
	gcvt(((y1+y2)/2+(y2-y1)/4),6,String);
	outtextxy((X1+X2)/2+8,(Y1+Y2)/2+(Y2-Y1)/4-2,String);
}

//实坐标中实数x转换到窗口坐标的值
int xToX(int X1,int X2,float x1,float x2,float x)
{
	return((int)(X1+(X2-X1)*(x-x1)/(x2-x1)));
}

//实坐标中实数y转换到窗口坐标的值
int yToY(int Y1,int Y2,float y1,float y2,float y)
{
	return((int)(Y1+(Y2-Y1)*(y-y1)/(y2-y1)));
}

//画插值函数曲线
void DrawFunctionCurve(int X1,int Y1,int X2,int Y2,
					   float x1,float y1,float x2,float y2,
					   POINT Table[],int n)
{
	int i;
	float Real_x_step;

	if(X2<=X1) {printf("error screen window!");return;}

	Real_x_step=(x2-x1)/(X2-X1);
	for(i=0;i<X2-X1;++i)
		line(xToX(X1,X2,x1,x2,x1+i*Real_x_step),
			 yToY(Y1,Y2,y1,y2,Lagrange(Table,n,x1+i*Real_x_step)),
			 xToX(X1,X2,x1,x2,x1+(i+1)*Real_x_step),
			 yToY(Y1,Y2,y1,y2,Lagrange(Table,n,x1+(i+1)*Real_x_step)),
			 3);
}

//在屏幕上显示一个看起来大一点的点
void PutBigDot(int x,int y,int color)
{
	setcolor(color);
	circle(x,y,2);
}

//画函数表中各点到屏幕上
void DrawPointsInTable(int X1,int Y1,int X2,int Y2,
					   float x1,float y1,float x2,float y2,
					   POINT Table[],int n)
{
	int i;

	for(i=0;i<n;++i)
		PutBigDot(xToX(X1,X2,x1,x2,Table[i].x),yToY(Y1,Y2,y1,y2,Table[i].y),14);
}

//画一点在屏幕上
void DrawPoint_xy(int X1,int Y1,int X2,int Y2,
				  float x1,float y1,float x2,float y2,
				  float x,float y)
{
	char Digits[60];
	int X,Y,tmp;

	X=xToX(X1,X2,x1,x2,x);
	Y=yToY(Y1,Y2,y1,y2,y);
	PutBigDot(X,Y,13);

	Digits[0]='(';
	gcvt(x,8,&Digits[1]);
	tmp=strlen(Digits);
	Digits[tmp]=',';
	Digits[tmp+1]=0;
	gcvt(y,8,&Digits[strlen(Digits)]);
	tmp=strlen(Digits);
	Digits[tmp]=')';
	Digits[tmp+1]=0;
	outtextxy(X+3,Y-3,Digits,13);
}

//画插值曲线和结果点
void DrawCurve_And_ShowResult(POINT Table[],int n,float x,float y)
{
	float Max_x,Max_y;
	float x1,y1,x2,y2;
	int X1,Y1,X2,Y2;

	Table[n].x=x;
	Table[n].y=y;
	Max_x=GetMaximum_x(Table,n+1);
	Max_y=GetMaximum_y(Table,n+1);

	x1=-Max_x*2;
	y1= Max_y*2;
	x2= Max_x*2;
	y2=-Max_y*2;
	X1=10;
	Y1=60;
	X2=630;
	Y2=450;

	DrawCoordinate(X1,Y1,X2,Y2,x1,y1,x2,y2);
	DrawFunctionCurve(X1,Y1,X2,Y2,x1,y1,x2,y2,Table,n);
	DrawPointsInTable(X1,Y1,X2,Y2,x1,y1,x2,y2,Table,n);
	DrawPoint_xy(X1,Y1,X2,Y2,x1,y1,x2,y2,x,y);
}

//函数表中是不是有x值相同的点,如有相同,则有错
int Check_x_is_same_as_former(POINT Table[],int i)
{
	int k;

	for(k=0;k<i;++k)
		if(fabs(Table[k].x-Table[i].x)<0.000001) return 1;
	return 0;
}

//输入函数表
void InputTable(POINT Table[],int *n)
{
	int i;

	printf("\nNumber of Points n=");
	scanf("%d",n);

	for(i=0;i<*n;++i)
	{
		printf("Input %dth point(x%d,y%d): ",i,i,i);
		scanf("%f,%f",&Table[i].x,&Table[i].y);
		if(Check_x_is_same_as_former(Table,i))
		{
			printf("\nThis x is the same as the former, ReInput!\n",i);
			--i;
		}
	}
}

//插值函数,这是核心函数,
float Lagrange(POINT Table[],int n,float x)
{
	float Result,tmp;
	int i,j;

	Result=0.0;
	for(i=0;i<n;++i)
	{
		tmp=1.0;
		for(j=0;j<n;++j)
		{
			if(j==i)continue;
			tmp*=((x-Table[j].x)/(Table[i].x-Table[j].x));
		}
		Result+=(tmp*Table[i].y);
	}
	return(Result);
}

void main()
{
	int n;
	float x,y;
	POINT Table[1000];

	//输入函数表
	InputTable(Table,&n);
	
	//输入一个x点,即要进行插值计算的点
	printf("\nInput x=");
	scanf("%f",&x);

	//计算Lagrange插值
	y = Lagrange(Table,n,x);

	//结果输出
	printf("\nL(%f)=%f",x,y);

	printf("\n\nPress any key to show Table,draw L(x) and point (x,y)\n");
	getch();

	//下面是有关画曲线和结果的图型化显示
	InitGraph();

	ShowTable(Table,n);
	DrawCurve_And_ShowResult(Table,n,x,y);

	outtextxy(1,468,"Press any key to quit!");
	getch();
	closegraph();
}

//注意
//用TC,BC作图形化程序时,要在菜单中把link库中的图型库选项打上X,否则有link时会报错

/*
	运行实例：(注意：当前目录下必须有egavga.bgi文件)

Number of Points n=3
Input 0th point(x0,y0): 100,10
Input 1th point(x1,y1): 121,11
Input 2th point(x2,y2): 144,12

Input x=115

L(115.000000)=10.722755

Press any key to show Table,draw L(x) and point (x,y)
*/