【C/C++】拉格朗日插值、牛顿插值、分段线性插值作业

本文只是完成作业，重点是比较三种插值过程和结果的特点（略），主函数（略），不考虑太多健壮性，也不批量插值。

创建项目

1. 使用VS2019创建C++控制台应用

   C++菜鸟教程

拉格朗日插值

公式

n次的拉格朗日插值多项式：$L_n(x)=\sum _{k=0}^n[y_k\prod _{i=0,i\ne k}^n(\frac{x-x_i}{x_k-x_i})]$

其中，已知点有n+1个，其x值依次为$x_0$、$x_1$、…$x_i$（或$x_k$）、…$x_n$，其y值同理。

当取n=1时，为线性插值；当取n=2时，为抛物插值。

思路

1.编程前，令上述公式中的n取n-1，于是公式为：$L_{n-1}(x)=\sum _{k=0}^{n-1}[y_k\prod _{i=0,i\ne k}^{n-1}(\frac{x-x_i}{x_k-x_i})]$

2.这样编写的程序函数中，接收的形参包含四个部分：

• x值的数组
• y值得数组
• 点数n
• 插值位置x

3.返回参数：

• 插值结果y

4.数据类型选择：double

5.为了对比不同插值方法的效率和结果，不建立额外的“缓存”

代码

参考博客1，可能有用的博客2（python版）

double lagrange(double arrX[], double arrY[], int n, double x)
{
	int k, i;
	double temp;
	double y = 0;
	for (k = 0; k < n; k++)
	{
		temp = 1;
		for (i = 0; i < n; i++)
		{
			if (i != k) {
				temp *= ((x - arrX[i]) / (arrX[k] - arrX[i]));
			}
		}
		y += arrY[k] * temp;
	}
	return y;
}

特点

优点：直观、简单、应用广泛。
缺点：当插值精度不够，增加新节点时，必须从头计算，不能利用已有结果。

牛顿插值

为克服拉格朗日插值的缺点，实现灵活增加插值节点，以节省运算次数。

参考博客1

公式

$N_n(x)=f(x_0)+\sum _{i=1}^n[\prod _{k=0}^{i-1}(x-x_i)]f[x_0,x_1,...,x_i]$

思路

1. 根据已知数据的个数n，建立差商表，表内共需要存储$\frac{n(n-1)}{2}$个差商值。
2. 根据差商表，将插值数据代入公式。

代码

计算差商表：

//计算差商表，n表示有n个节点,n>1，f为差商表数组
void calcChaShang(double x[], double y[], double f[], int n)
{ 
	int i, j, k = 0;
	//一阶差商计算
	for (i = 0; i < n - 1; i++)
	{
		f[k++] = (y[i] - y[i + 1]) / (x[i] - x[i + 1]);
	}

	//在一阶差商的基础上，计算到n-1阶
	for (i = 1; i < n - 1; i++)
	{
		//j用来遍历第i阶的所有差商(共有n-i个)
		for (j = 0; j < n - i - 1; j++)
		{
			//调试时所用
			/*double y1 = f[k - (n - i)];
			double y2 = f[k - (n - i) + 1];
			double x1 = x[j];
			double x2 = x[j + i + 1];*/
			//由于k自增会影响y[]，所以y[]中的j可以去掉，否则错误写法：
			//f[k] = (f[k - (n - i) + j] - f[k - (n - i) + j + 1]) / (x[j] - x[j + i + 1]);
			f[k] = (f[k - (n - i)] - f[k - (n - i) + 1]) / (x[j] - x[j + i + 1]);
			k++;
		}
	}
}

基于差商表的牛顿插值：

double newvalue(double* xArr, double* yArr, double* f, int n, double x)
{ 
	//将连乘结果存储到数组中（multiplication）
	double *m = (double*)malloc(sizeof(double) * n);
	m[0] = 1.0;
	int i = 0;
	for (i = 0; i < n - 1; i++)
	{
		m[i + 1] = m[i] * (x - xArr[i]);
	}
		
	//最终计算
	int k = 0;//通过此下标，访问差商数组
	double y = yArr[0];
	for (i = 1; i < n; i++) 
	{
		y += m[i] * f[k];
		k += n - i;
	}
	return y;
}

线性分段插值

公式

在区间[a,b]，划分$a=x_0<x_1<x_2<...<x_n=b$，对于[a,b]之间的任一小区间$[x_{j-1},x_j]$，在该小区间上作线性插值：

$$f(x)\approx L_1(x)=y_{j-1}\frac{x-x_j}{x_{j-1}-x_j}+y_j\frac{x-x_{j-1}}{x_j-x_{j-1}}$$

思路

1. 已知数据的x值需要从小到大排列。
2. 插值前先找到被插数据x值所在的分段，通过比较$x$和$x_{j-1}$的值，找到其所在的小区间$[x_{j-1},x_j]$中的j的值。
3. 找到j，代入公式。

代码

double linear(double xArr[], double yArr[], int n, double x) 
{
	int j;
	double y = 0;
	bool noFound = true;

	//如果是左边外插
	if (x <= xArr[0])
	{
		j = 0;
	}
	//如果是右边外插
	else if (x >= xArr[n - 1])
	{
		j = n - 2;
	}
	//内插
	else
	{
		for (j = 1; j < n; j++)
		{
			//找到所在的分段
			if (x <= xArr[j])
			{
				j--;
				break;
			}
		}
	}

	y = yArr[j] * ((x - xArr[j + 1]) / (xArr[j] - xArr[j + 1]))
		+ yArr[j + 1] * ((x - xArr[j]) / (xArr[j + 1] - xArr[j]));
	return y;
}

参考

公式来自《数值计算方法及其应用》（朱长青编著，科学出版社）