武汉理工大学-数值分析-(1)插值方法

最新推荐文章于 2024-06-28 16:20:27 发布
最新推荐文章于 2024-06-28 16:20:27 发布
阅读量3.8k 收藏 37
点赞数 2
分类专栏： 数值分析 文章标签： c++ 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_43381409/article/details/104396075
版权

文章目录

实验目标

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

编程语言与扩展库

语言：C++语言
导入扩展库：Eigen库
在这里插入图片描述
在项目名称处 右键->属性->VC++目录->包含目录->选择eigen模块所在文件夹路径，即可导入模块；
eigen资源文件可以在eigen官网上进行下载：http://eigen.tuxfamily.org/index.php?title=Main_Page

拉格朗日插值法

//拉格朗日插值法
//参数：插值点 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;
}

计算公式：对于 拟合曲线 y = a + bx
在这里插入图片描述

源代码整合

#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文件 中，
通过 函数 ReadPointData 读取点数据。

points1.txt：（拉格朗日插值）
在这里插入图片描述
points2.txt：（牛顿插值）
在这里插入图片描述
points3.txt：（最小二乘线性拟合）
在这里插入图片描述
结果：
在这里插入图片描述

写在最后

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

-Kingzy-
关注
专栏目录
四元数插值Eigen源码解析
Volador的博客
05-10 1018
【代码】四元数插值Eigen源码解析。
利用Eigen进行B-Spline 曲线插值运算
小蚊子的博客
01-12 6146
B-Spline curve即B样条曲线，1972 年Gordon、Riesenfeld 等人提出，为了解决Bezier曲线的不足（Bezier 曲线不能做局部修改；Bezier 曲线的拼接比较复杂）。B样条曲线的计算原理可参考这里，或者看中文的更详细的介绍可参考《计算机辅助几何设计与非均匀有理B样条》这本书。     Eigen库unsupported部分里有实现B-Spline 曲线插值运算
Eigen中关于四元数的常用操作
最新发布
点云兔子的博客
06-28 653
Eigen库中的四元数类提供了强大的工具来处理三维旋转和插值。通过Eigen库，我们可以轻松实现四元数的创建、操作和插值，从而简化复杂的三维旋转问题。
武汉理工大学数值分析试题
06-20
这是最近几年武汉理工大学数值分析的一次考试的考试试题，
武汉理工大学-数值分析-(2)数值积分与微分
-Kingzy- 的博客
02-19 1888
文章目录实验目标编程语言与扩展库复化梯形积分龙贝格算法（Romberg）源代码整合运行结果写在最后 实验目标 使用编程语言实现以下算法： 1.复化梯形公式 的 自动控制误差 算法； 2.龙贝格积分 算法. 计算某个区间 [a,b] 上的积分： 编程语言与扩展库 语言：C++语言 导入扩展库：Eigen库 在项目名称处 右键->属性->VC++目录->包含目录->选择ei...
武汉理工大学数值分析欧拉
05-17
功能好，C++编码
武汉理工大学-数值分析-(3)线性代数方程组的数值解法
-Kingzy- 的博客
02-19 1509
文章目录实验目标编程语言与扩展库高斯消元法（Gauss）列主元高斯消元法（Gauss）列主元LU直接分解法Jacobi迭代法Gauss-Seidel迭代法写在最后 实验目标 用编程语言编程实现以下算法： 1.用 高斯(Gauss)消元法 求n阶线性方程组的解。 2.用 列主元高斯(Gauss)消元法 求n阶线性方程组的解。 3.用 列主元LU直接分解法 求n阶线性方程组的解。 4.用 Jacobi...
数值分析复习提纲-武汉理工大学版
05-22
数值分析复习提纲-武汉理工大学版带目录更正版，纯手打
武汉理工大学数值分析历年真题与答案解析，陈晓江数值分析课后习题答案
05-28
本资源集合了武汉理工大学数值分析课程的相关学习材料，包括历年真题及答案解析，以及陈晓江教授编写的《数值分析》教材的课后习题解答，对于学习者来说是一份宝贵的参考资料。 首先，"数值分析考试真题.pdf"可能...
武汉理工大学数值分析实验代码
10-08
在压缩包中的“武汉理工大学数值分析实验报告及代码”文件，我们可以深入研究每一个实验，分析其算法的效率和精度，同时对比不同数值方法的优缺点。这将有助于深化对数值分析理论的理解，并提高编程技能。通过这样的...
mathtoolbox：用Eigen用C ++ 11编写的数学工具（插值，降维，优化等）
01-31
数学工具箱 用C ++ 11和编写的数学工具（插值，降维，优化等）。 演算法 分散数据插值和函数逼近 降维和低维嵌入 数值优化 线性代数 实用工具 依存关系 主图书馆 本征 brew install eigen / sudo apt install libeigen3-dev ） Python绑定 pybind11 （包含在gitsubmodule中） 例子 优化测试功能 （包含在gitsubmodule中） 计时器 （包含在gitsubmodule中） 用作C ++库 mathtoolbox使用CMake 来构建源代码。 例如，可以通过以下方式构建该库 git clone https://github.com/yuki-koyama/mathtoolbox.git --recursive cd mathtoolbox mkdir build cd build cmake ../ make 并且可以选择通过以下方式将其安装到系统中 make install 当CMake参数MATHTOOLBOX_BUILD_EXAMPLES设置为ON ，还将构建示例应用程序。 （默认设置为O
用Eigen用C ++ 11编写的数学工具（插值，降维，优化等）-C/C++开发
05-27
用Eigen mathtoolbox用C ++ 11编写的数学工具（插值，降维，优化等）用C ++ 11和Eigen编写的数学工具（插值，降维，优化等）。 算法散点数据插值和函数逼近rbf插值：径向基函数（RBF）插值高斯过程回归：高斯过程回归（GPR）降维和低维嵌入classic-mds：经典多维缩放（MDS）som ：自组织映射（SOM）数值优化
数值分析(实验一、二).zip
12-17
武汉理工大学数值分析实验一（插值方法：拉格朗日插值多项式、牛顿插值、线性函数拟合）代码和实验报告、实验二（数值积分与微分与常微分方程初值问题的数值解法matlab：复化梯形积分、Romberg积分、改进的欧拉...
计算方法 | 给定插值点 (xi , yi), i = 0, 1, 2, · · · , n，可分别构造 Lagrange 插值多项式和 Newton 插值多项式，证明两者相同并说明各自具有的特点？
黄佳俊的博客
12-06 1279
证明两者相同： 各自特点：
Eigen
lne_fly的博客
11-25 136
http://eigen.tuxfamily.org/index.php?title=Main_Page#Documentation
C++线性代数运算库eigen3的使用，空间欧式变换的表示和运算，旋转四元数的球面插值
qq_41253960的博客
03-25 5378
c++库eigen, 矩阵运算
Eigen线性代数库学习大全
几度春风里的博客
11-17 2761
Eigen是C++的线性代数库，能提供有关矩阵的线性代数运算，还包含解方程等功能。
直线段-两点线性插值及可视化
baidu_38227555的博客
04-10 1002
（博客园最近不知道怎么了，把自己的随笔搬一些过来CSDN这边） 已知直线段AB，线性插值求插值点C的公式为：C =A*（1.0 -t） +B*t，，其中t的范围为 [0.0, 1.0], 0值时候为A点，1值时候为B点 下面是用QPainter可视化 GrphWidget.h #ifndef GRPHWIDGET_H #define GRPHWIDGET_H #include <QWidget> #include "ui_grphwidget.h" class Grph...
数值分析：插值与逼近方法详解
"数值分析第六章主要探讨了插值与逼近的方法，包括Lagrange插值、Newton插值、Hermite插值、分段插值以及三次样条插值等多种技术，同时也提到了拟合方法如最小二乘法。这些方法在处理实验数据时起到重要作用，能够将...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值