【数值分析课程实验-01】插值方法

实验目标

使用编程语言实现以下算法：
1. 已知插值节点序列，使用 拉格朗日 (Lagrange) 插值 多项式计算函数在点的近似值；
2. 已知插值节点序列，使用 牛顿 (Newton) 插值 多项式计算函数在点的近似值；
3. 使用 线性函数 的 最小二乘法 对给定数据进行 拟合。

编程语言与扩展库

编程语言：C++ 语言
IDE 环境：Visual Studio 2017
扩展库：Eigen 库

在项目名称处 右键 >> 属性 >> VC++目录 >> 包含目录 >> 选择 Eigen 模块所在文件夹路径，导入模块；
Eigen 资源文件可以在 Eigen 官网上进行下载。

实现代码

拉格朗日插值法

// 拉格朗日插值法
// 参数：插值点 points，计算点 x，插值点个数 n
double Lagrange(vector<Vector2d> points, double x, int n)
{
	// 运算结果
	double result = 0;
	// 插值基函数
	double L;

	for (int k = 0; k < n; k++) {
		// 计算插值基函数
		L = 1;
		for (int j = 0; j < n; j++) {
			// 叠乘 (x-xj) / (xk-xj)
			if (j != k) {
				L = ((x - points[j][0]) / (points[k][0] - points[j][0]))*L;
			}
		}
		// 叠加 Lk * yk
		result += L * points[k][1];
	}
    
	return result;
}

计算公式：

$$p_k\left( x \right) =\prod_{i\in B_k}{\frac{x-x_i}{x_k-x_i}}\ ,\ \ \ \ L_n\left( x \right) =\sum_{j=0}^{n-1}{y_jp_j\left( x \right)}$$

算法流程图：

牛顿插值法

算法通过 计算差商表 实现：

// 牛顿插值法
double Newton(vector<Vector2d> points, double x, int n)
{
	// 运算结果
	double result = 0;

	// 差商表
	MatrixXd DQTable(n, n + 1);
	DQTable.setZero();

	// 构造差商表 (i 为行坐标，j 为列坐标)
	for (int i = 0; i < n; i++) {
		// 第 0 列为各插值点的 x 值
		DQTable(i, 0) = points[i][0];
	}
	for (int i = 0; i < n; i++) {
		// 第 1 列为各插值点的 y 值
		DQTable(i, 1) = points[i][1];
	}
	// 计算各阶差商
	for (int j = 2; j < n + 1; j++) {
		for (int i = j - 1; i < n; i++) {
			DQTable(i, j) = (DQTable(i, j - 1) - DQTable(i - 1, j - 1)) / (DQTable(i, 0) - DQTable(i - j + 1, 0));
		}
	}

	// w = (x-x0)(x-x1)...(x-xn-1) 叠乘实现
	double w = 1;
	for (int k = 0; k < n; k++) {
		if (k > 0) {
			w = w * (x - points[k - 1][0]);
		}
		// 叠加 w * k 阶差商 (差商表对角线)
		result += w * DQTable(k, k + 1);
	}

	return result;
}

计算公式：

$$f\left( x,x_0,x_1,...,x_n \right) =\frac{f\left( x_0,x_1,...,x_n \right) -f\left( x,x_0,...,x_{n-1} \right)}{x_n-x}$$ $$f\left( x \right) =f\left( x_0 \right) +\left( x-x_0 \right) f\left( x,x_0 \right) +\left( x-x_0 \right) \left( x-x_1 \right) f\left( x_0,x_1,x_2 \right)$$ $$+\ ...\ +\left( x-x_0 \right) \left( x-x_1 \right) ...\left( x-x_n \right) f\left( x_0,x_1,...,x_n \right)$$

最小二乘线性拟合

算法采用 法方程法 计算：

// 一次多项式曲线拟合  y = a + b*x
// 返回参数 a、b 组成的向量
Vector2d CurveFit(vector<Vector2d> points, int n)
{
	double sum_x = 0;    // x 总和
	double sum_x2 = 0;    // x^2 总和
	double sum_y = 0;    // y 总和
	double sum_xy = 0;    // x*y 总和

	Matrix2d A;    // 2*2 矩阵
	Vector2d b;    // 向量
	Vector2d x;    // 求解向量

	for (int i = 0; i < n; i++) {
		sum_x += points[i][0];
		sum_x2 += pow(points[i][0], 2);
		sum_y += points[i][1];
		sum_xy += points[i][0] * points[i][1];
	}

	A(0, 0) = n;
	A(1, 0) = A(0, 1) = sum_x;
	A(1, 1) = sum_x2;

	b(0) = sum_y;
	b(1) = sum_xy;

	// 求解方程 Ax = b
	x = A.lu().solve(b);

	return x;
}

计算公式：对于 拟合曲线 y=a+bx ，按照下列方程求解列向量 x 。

$$\left[ \begin{matrix} n& \sum{x}\\ \sum{x}& \sum{x^2}\\ \end{matrix} \right] \left[ \begin{array}{c} a\\ b\\ \end{array} \right] =\left[ \begin{array}{c} \sum{y}\\ \sum{xy}\\ \end{array} \right]$$ $$Ax\ =\ b$$

源代码整合

#include <Eigen/Dense>
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>

using namespace Eigen;
using namespace std;

// 拉格朗日插值法
// 参数：插值点 points，计算点 x，插值点个数 n
double Lagrange(vector<Vector2d> points, double x, int n)
{
	// 运算结果
	double result = 0;
	// 插值基函数
	double L;

	for (int k = 0; k < n; k++) {
		// 计算插值基函数
		L = 1;
		for (int j = 0; j < n; j++) {
			// 叠乘 (x-xj) / (xk-xj)
			if (j != k) {
				L = ((x - points[j][0]) / (points[k][0] - points[j][0]))*L;
			}
		}
		// 叠加 Lk * yk
		result += L * points[k][1];
	}

	return result;
}

// 牛顿插值法
double Newton(vector<Vector2d> points, double x, int n)
{
	// 运算结果
	double result = 0;

	// 差商表
	MatrixXd DQTable(n, n + 1);
	DQTable.setZero();

	// 构造差商表 (i 为行坐标，j 为列坐标)
	for (int i = 0; i < n; i++) {
		// 第 0 列为各插值点的 x 值
		DQTable(i, 0) = points[i][0];
	}
	for (int i = 0; i < n; i++) {
		// 第 1 列为各插值点的 y 值
		DQTable(i, 1) = points[i][1];
	}
	// 计算各阶差商
	for (int j = 2; j < n + 1; j++) {
		for (int i = j - 1; i < n; i++) {
			DQTable(i, j) = (DQTable(i, j - 1) - DQTable(i - 1, j - 1)) / (DQTable(i, 0) - DQTable(i - j + 1, 0));
		}
	}

	// w = (x-x0)(x-x1)...(x-xn-1) 叠乘实现
	double w = 1;
	for (int k = 0; k < n; k++) {
		if (k > 0) {
			w = w * (x - points[k - 1][0]);
		}
		// 叠加 w * k 阶差商 (差商表对角线)
		result += w * DQTable(k, k + 1);
	}

	return result;
}

// 一次多项式曲线拟合  y = a + b*x
// 返回参数 a、b 组成的向量
Vector2d CurveFit(vector<Vector2d> points, int n)
{
	double sum_x = 0;    // x 总和
	double sum_x2 = 0;    // x^2 总和
	double sum_y = 0;    // y 总和
	double sum_xy = 0;    // x*y 总和

	Matrix2d A;    // 2*2 矩阵
	Vector2d b;    // 向量
	Vector2d x;    // 求解向量

	for (int i = 0; i < n; i++) {
		sum_x += points[i][0];
		sum_x2 += pow(points[i][0], 2);
		sum_y += points[i][1];
		sum_xy += points[i][0] * points[i][1];
	}

	A(0, 0) = n;
	A(1, 0) = A(0, 1) = sum_x;
	A(1, 1) = sum_x2;

	b(0) = sum_y;
	b(1) = sum_xy;

	// 求解方程 Ax = b
	x = A.lu().solve(b);

	return x;
}

// 读取点数据文件的函数
// 参数：文件路径 path，计算点 x0
vector<Vector2d> ReadPointData(const string &path, double &x0)
{
	// 点集
	vector<Vector2d> points;

	// 建立读入流对象
	ifstream ifs(path);
	if (!ifs.is_open()) return points;

	// 先读入计算点 x 坐标
	ifs >> x0;

	// 循环读入点
	double x, y;
	while (!ifs.eof()) {
		ifs >> x >> y;
		points.push_back(Vector2d(x, y));
	}

	return points;
}

int main()
{
	// 插值点集
	vector<Vector2d> points;   
	// 曲线拟合参数 a、b 
	Vector2d v;    
	// 运算点 x 坐标
	double x;       

	// 拉格朗日插值
	points = ReadPointData("F://points1.txt", x);
	cout << endl << "拉格朗日插值：" << "f(" << x << ") = " << Lagrange(points, x, points.size()) << endl;

	// 牛顿插值
	points = ReadPointData("F://points2.txt", x);
	cout << endl << "牛顿插值：" << "f(" << x << ") = " << Newton(points, x, points.size()) << endl;

	// 曲线拟合
	points = ReadPointData("F://points3.txt", x);  v = CurveFit(points, points.size());
	cout << endl << "曲线拟合：" << "y = " << v[0] << " + " << v[1] << " * x" << endl;

	return 0;
}

运行结果

三个算法的实验数据分别存储在 F 盘的三个 txt 文档中，
分别为：points1.txt (拉格朗日插值)、points2.txt (牛顿插值)、points3.txt (最小二乘线性拟合)。

运行结果：

写在最后

声明：本专栏内容来源于武汉理工大学 2019-2020 学年数值分析方法课程实验，仅供学习参考。
如有错误或不足地方，还请指出，祝愿读者能够在编程之路上不断进步！