武汉理工大学-数值分析-(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线性代数库学习大全
几度春风里的博客
11-17 2763
Eigen是C++的线性代数库，能提供有关矩阵的线性代数运算，还包含解方程等功能。
用Eigen用C ++ 11编写的数学工具（插值，降维，优化等）-C/C++开发
05-27
用Eigen mathtoolbox用C ++ 11编写的数学工具（插值，降维，优化等）用C ++ 11和Eigen编写的数学工具（插值，降维，优化等）。 算法散点数据插值和函数逼近rbf插值：径向基函数（RBF）插值高斯过程回归：高斯过程回归（GPR）降维和低维嵌入classic-mds：经典多维缩放（MDS）som ：自组织映射（SOM）数值优化
四元数插值Eigen源码解析
Volador的博客
05-10 1020
【代码】四元数插值Eigen源码解析。
利用Eigen进行B-Spline 曲线插值运算
小蚊子的博客
01-12 6147
B-Spline curve即B样条曲线，1972 年Gordon、Riesenfeld 等人提出，为了解决Bezier曲线的不足（Bezier 曲线不能做局部修改；Bezier 曲线的拼接比较复杂）。B样条曲线的计算原理可参考这里，或者看中文的更详细的介绍可参考《计算机辅助几何设计与非均匀有理B样条》这本书。     Eigen库unsupported部分里有实现B-Spline 曲线插值运算
Eigen中关于四元数的常用操作
最新发布
点云兔子的博客
06-28 660
Eigen库中的四元数类提供了强大的工具来处理三维旋转和插值。通过Eigen库，我们可以轻松实现四元数的创建、操作和插值，从而简化复杂的三维旋转问题。
C++线性代数运算库eigen3的使用，空间欧式变换的表示和运算，旋转四元数的球面插值
qq_41253960的博客
03-25 5378
c++库eigen, 矩阵运算
使用Eigen函数库做球面四边形插值的方法
zgce的专栏
04-26 3002
正在将一段D3D的代码改造成OpenGL的代码，在原来D3D中使用了D3DXQuaternionSquadSetup与D3DXQuaternionSquad两个D3D的数学函数。 其中D3DXQuaternionSquadSetup这个函数的MSDN解释如下： void D3DXQuaternionSquadSetup( __in D3DXQUATERNION *pAOut,
手撕自动驾驶算法——四元素插值
xiaohu的博客
09-23 2258
四元数插值常见的有线性插值、球面线性插值等 线性插值（Linear Interpolation，简称Lerp） 四元数表示旋转时是单位四元数，这种插值方式，相当于我们是沿着一条直线（也就是圆上的一个弦）进行插值的，这样插值出来的四元数不是单位四元数 归一化线性插值（Normalized LinearInterpolation，简称Nlerp) 有其他问题的。如下图所示，在同等时间内， vt ...
武汉理工大学-数值分析-(3)线性代数方程组的数值解法
01-06
1.用 高斯(Gauss)消元法 求n阶线性方程组的解。 2.用 列主元高斯(Gauss)消元法 求n阶线性方程组的解。 3.用 列主元LU直接分解法 求n阶线性方程组的解。 4.用 Jacobi迭代法 求n阶线性方程组的解。 5.用 Gauss-Seidel...
武汉理工大学数值分析历年真题与答案解析，陈晓江数值分析课后习题答案
05-28
本资源集合了武汉理工大学数值分析课程的相关学习材料，包括历年真题及答案解析，以及陈晓江教授编写的《数值分析》教材的课后习题解答，对于学习者来说是一份宝贵的参考资料。 首先，"数值分析考试真题.pdf"可能...
数值分析复习提纲-武汉理工大学版
05-22
这个复习提纲来自武汉理工大学，旨在帮助学生系统地理解和掌握数值分析的基本概念、理论和方法。下面我们将深入探讨其中涉及的主要知识点。 1. **数值线性代数** - **矩阵运算**：包括矩阵的加减、乘法和求逆，...
武汉理工大学数值分析实验代码
10-08
在压缩包中的“武汉理工大学数值分析实验报告及代码”文件，我们可以深入研究每一个实验，分析其算法的效率和精度，同时对比不同数值方法的优缺点。这将有助于深化对数值分析理论的理解，并提高编程技能。通过这样的...
武汉理工大学数值分析拉格朗日插值C++代码
05-17
### 武汉理工大学数值分析拉格朗日插值C++代码解析 #### 概述 在数学和计算机科学领域，插值是一种通过已知数据点来估算未知数据点的方法。拉格朗日插值法是一种常用的多项式插值方法，它能够通过一系列离散的数据...
Eigen矩阵运算库使用
m0_37870385的博客
05-22 2496
文章目录lpl_plp​ norm／p范数计算colwise()/rowwise()array()blocklpl_plp​ norm／p范数计算 L2范数 squareNorm()，等价于计算vector的自身点积，norm()返回squareNorm的开方根。 这些操作应用于matrix，norm() 会返回Frobenius或Hilbert-Schmidt范数。 如果你想使用其他Lp范数，可以使用lpNorm< p >()方法。p可以取Infinity，表示L∞范数。 int main()
对3d坐标和四元数做线性插值计算c++实现
栾金鹿
06-06 1753
对3d坐标和四元数做线性插值计算，c++实现
武汉理工大学-数值分析-(4)方程求根的数值方法
-Kingzy- 的博客
02-19 1255
文章目录实验目标编程语言与扩展库二分法求根牛顿法求根弦截法求根写在最后 实验目标 用编程语言编程实现以下算法： 1、用 二分法 求解方程的根； 2、用 牛顿(Newton)迭代法 求解方程的根； 3、用 弦截法 求解方程的根. 编程语言与扩展库 语言：Python 扩展模块：numpy 二分法求根 # 函数 f(x) def f(x): y = x**3 - x - 1 re...
直线段-两点线性插值及可视化
baidu_38227555的博客
04-10 1003
（博客园最近不知道怎么了，把自己的随笔搬一些过来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...
三次样条插值c2连续 cubic C2 spline
home180490的博客
12-17 1118
cubic C2 spline 主要工作其实就是求解a,b,c,d四个参数 通过以上2次可微等公式，可解得参数为： #include "cubicSpline.h" namespace cubicSpline { std::vector<double> vec_diff(std::vector<double> input) { std::vector<double> output; for (unsigned int i = 1; i < inpu
数值分析：插值与逼近方法详解
"数值分析第六章主要探讨了插值与逼近的方法，包括Lagrange插值、Newton插值、Hermite插值、分段插值以及三次样条插值等多种技术，同时也提到了拟合方法如最小二乘法。这些方法在处理实验数据时起到重要作用，能够将...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值