数值分析实验相关

综合说明：

在写这次实验的时候，程序的图形方式选用的是SDL图形库，关于SDL图形库的使用，可以网上查询相关资料，对于一般的2D图形显示还是很简单好上手的。

这些实验的代码我都发在了GitHub：https://github.com/FieldSoft-HelloClyde/NumericalAnalysisExperiment/

也遇到一些问题，我也分享一下解决办法：

1.配置完的SDL项目显示控制台：

需要在项目属性-配置属性-链接器-系统 选项卡中把子系统设置为控制台。

2.VS2012程序不能再XP系统上运行：

首先需要对VS2012打上最新的补丁。访问https://www.visualstudio.com/downloads/download-visual-studio-vs#d-additional-software，选择VS2012更新X（目前是5）.

下载安装后项目属性-配置属性-常规选项卡下中的平台工具集可以选择Windows XP相关选项

但是这样往往还不够，可能在XP运行的时侯还会出现缺少msvcp110.dll的问题。

这就需要安装VS2012的运行库

http://www.microsoft.com/zh-cn/download/details.aspx?id=30679

下载安装即可。

3.调试时候遇到程序相对路径不正确的问题

解决方法是在项目属性-配置属性-调试-工作目录 选择生成的程序根目录

4.SDL 鼠标滚轮使用

这里贴一段代码大概就明白了。

while(SDL_PollEvent(&event))
		{
			if (event.type == SDL_QUIT){
				While_Quit = 4;
			}
			else if (event.type == SDL_MOUSEMOTION){
				if (MoveState){
					int dx,dy;
					dx = event.motion.x - Mouse_Old_X;
					dy = event.motion.y - Mouse_Old_Y;
					X_pos += dx;
					Y_pos += dy;
					Mouse_Old_X = event.motion.x;
					Mouse_Old_Y = event.motion.y;
				}
			}
			else if (event.type == SDL_MOUSEBUTTONDOWN){
				if (event.button.button == SDL_BUTTON_LEFT){
					if (event.motion.x >= 20 && event.button.x <= 65 && event.button.y >= 80 && event.button.y <= 100){
						X_pos = Windows_Width / 2;
						Y_pos = Windows_Height / 2;
						PixRate = 1.0;
					}
					else{
						MoveState = true;
						Mouse_Old_X = event.motion.x;
						Mouse_Old_Y = event.motion.y;
					}
				}
			}
			else if (event.type == SDL_MOUSEBUTTONUP){
				if (event.button.button == SDL_BUTTON_LEFT){
					MoveState = false;
				}
			}
			else if (event.type == SDL_MOUSEWHEEL){
				if (event.wheel.y > 0){
					//放大
					PixRate *= 2;
				}
				else if (event.wheel.y < 0){
					PixRate /= 2;
				}
			}
         }

当y值大于0时表示滚轮向上滚动，小于0表示滚轮向下滚动。

1.Lagrange插值多项式

公式：

关键算法：

double Lagrange_n(double xx){
	double x[] = {0,10,20,30,40,50,60,70,80,90,100,110,120};
	double y[] = {5,1,7.5,3,4.5,8.8,15.5,6.5,-5,-10,-2,4.5,7};
	return Lagrange(13,x,y,xx);
}

double Lagrange(int n,double *x,double *y,double xx){
	double addResult = 0;
	int k = 0;
	for (k = 0;k < n;k ++){
		double mulResult = 1;
		int j = 0;
		for (j = 0;j < n;j ++){
			if (j == k) continue;
			mulResult *= (xx - x[j]) / (x[k] - x[j]);
		}
		addResult += y[k] * mulResult;
	}
	return addResult;
}

其中数组x,y就是测试数据

运算结果：

所画函数就是Ln(x)，红点就是测试数据。可能看不清。

2.Newton插值多项式

公式：

其中f[x0,x1,....xn]为差商公式

具体公式为：（来自百度百科）

差商（difference quotient)的定义：

函数f(x)在两个互异点xi,xj处的一阶差商定义为：

f[xi,xj]=[f(xi)-f(xj)]/(xi-xj)=(yi-yj)/(xi-xj)

一阶差商的差商称为函数f(x)在xi,xj,xk(i,j,k互异)点的二阶差商：

f[xi,xj,xk]={f[xi,xj]-f[xj,xk]}/(xi-xk)

一般的，k-1阶差商的差商称为函数f(x)在n+1个互异点x0,x1,...,xk上的k阶差商：

f[x0,x1,...,xk]={f[x0,x1,...,xk-1]-f[x1,x2,...,xk]}/(x0-xk)

特别的，规定零阶差商f[xk]=f(xk)=yk

我们使用递归来定义它。

相关算法：

double Newton_n(double xx){
	double x[] = {0,10,20,30,40,50,60,70,80,90,100,110,120};
	double y[] = {5,1,7.5,3,4.5,8.8,15.5,6.5,-5,-10,-2,4.5,7};
	for (int i = 0;i < 13;i ++){
		for (int j = 0;j < 13;j ++){
			DiffQuoTable[i][j] = 0;
			DiffQuoExist[i][j] = false;
		}
	}
	return Newton(13,x,y,xx);
}

/*
原本效率较低
优化001：
	数组储存临时内容，减少重复计算
*/
double DiffQuotient(double *x,double *y,int xstart,int xx){
	if (DiffQuoExist[xstart][xx]){
		return DiffQuoTable[xstart][xx];
	}
	else{
		if (xx == 1){
			return (y[xstart] - y[xstart + 1]) / (x[xstart] - x[xstart + 1]);
		}
		else if (xx == 0){
			return y[xstart + xx];
		}
		else{
			DiffQuoTable[xstart][xx] = (DiffQuotient(x,y,xstart,xx - 1) - DiffQuotient(x,y,xstart + 1,xx - 1)) / (x[xstart] - x[xstart + xx]);
			DiffQuoExist[xstart][xx] = true;
			return DiffQuoTable[xstart][xx];
		}
	}
}

double Newton(int n,double *x,double *y,double xx){
	double result = y[0];
	int index = 0;
	for (index = 0;index < n - 1;index ++){
		double anresult;
		anresult = DiffQuotient(x,y,0,index + 1);
		int index2 = 0;
		for (index2 = 0;index2 <= index;index2 ++){
			anresult *= (xx - x[index2]);
		}
		result += anresult;
	}
	return result;
}

对于差商算法，只使用递归会导致效率极低。经过分析发现，差商计算过程中有很多中间步骤是重复计算的，于是使用数组来临时储存计算结果，来对其优化。效果很不错。

计算结果：

3.Hermite插值多项式

公式：

算法：

double Hermite_n(double xx){
	double x[] = {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0};
	double y[] = {0.904837,0.818731,0.740818,0.670320,0.606531,
				0.548812,0.496585,0.449329,0.406570,0.367879};
	double m[] = {-0.904837,-0.818731,-0.740818,-0.670320,-0.606531,
				-0.548812,-0.496585,-0.449329,-0.406570,-0.367879};
	return Hermite(10,x,y,m,xx);
}

double Hermite(int n,double *x,double *y,double *m,double xx){
	double Result = 0;
	for (int j = 0;j < n;j ++){
		double lj = 1;
		for (int i = 0;i < n;i ++){
			if (i == j) continue;
			lj *= (xx - x[i]) / (x[j] - x[i]);
		}
		double ljs = 0;
		for (int k = 0;k < n;k ++){
			if (k == j) continue;
			ljs += (double)1.0 / (x[j] - x[k]);
		}
		double aj;
		aj = (1 - 2 * (xx - x[j]) * ljs) * lj * lj;
		double bj;
		bj = (xx - x[j]) * lj * lj;
		Result += y[j] * aj + m[j] * bj;
	}
	return Result;
}

运算结果：