数值计算之函数逼近

目录

一、实验内容

二、实验过程

(一)概念的准备

1.函数逼近问题

2.插值法/内插法

3.基本思路

(二)拉格朗日插值多项式

1.基函数的构造

2.拉格朗日插值多项式的构造

3.代码的实现

 4.特别的函数

 5.优缺点

(三)牛顿插值

0.前言

1.插值基函数的优化

 2.差商----各项系数

3.代码的实现

4.优缺点

(四) 埃尔米特插值

1.简介

 2.基函数的确定

3.代码

4.优缺点

(五)西楚霸王的自编之范德蒙德待定系数法

1.原理

 2.代码

3.优缺点

三、实验结果

1.主函数

2.运行结果

3.四种方法的比较

四、实验感受


一、实验内容

(1) 给定函数e^-x在x=1,2,3时的值入下表。编写用拉格朗日法和牛顿法,分别用1次和2次插值计算e^-2.6的值。

(2) 查阅文献整理插值的其它方法,并至少编程实现其中的一种算法。结合(1)的结果对比分析几种算法的性能。

二、实验过程

(一)概念的准备

1.函数逼近问题

当函数只在有限点集上给定函数值,要在包含该点集的区间上用公式给出函数的简单表达式.或者是在区间上用简单函数逼近已知复杂函数的问题,就是函数逼近问题。

2.插值法/内插法

构造某个函数作为不便于处理或计算的函数的近似,然后通过处理简单函数获得近似结果,当要求近似函数取给定的离散数据时,这种处理方法称为插值法。

插值法是使用插值多项式来逼近函数,唯一的要求是插值函数和被插函数在插值节点上的函数值相同。

如果插值函数是次数不超过n的代数多项式(如上图),就称为多项式差值。

如果插值函数是分段多项式,就称为分段差值。

如果插值函数是三角多项式,就称为三角差值。

3.基本思路

构造一个比较简单的函数使该函数可以通过所以已知节点。

(二)拉格朗日插值多项式

1.基函数的构造

n=1时,求一次多项式:y=ax+b(代入x0,y0和x1,y1)

化简消去a、b,得到基函数l0、l1(这是小写的L);

观察两个基函数,有趣的是他们类似于开关的作用,只有在对应的xi处,基函数的值才为1,否则分子为0(开关关闭)。

推广到n个插值函数,因而拉格朗日插值法构造了n个插值基函数,第i个基函数保证在xi上为1,在其他xj上均为零。然后再将这些基函数分别乘以yi相加,所得函数便能满足要求。

2.拉格朗日插值多项式的构造

所构造的基函数 l 和函数 L :

3.代码的实现

//拉格朗日插值法,index表示插值点
double Lagrange(int n, double x[], double y[], double index)
{
	double res = 0;//记录y值 
	for (int i = 0; i < n; i++)//n
	{
		double temp = 1;//记录每一个小l 
		for (int j = 0; j < n; j++)
		{
			if (i != j)//n-1
				temp *= (index - x[j]) / (x[i] - x[j]);//这个是循环计算小l的 
		}
		res += temp * y[i];//累加每一项系数和基函数的积 
	}
	return res;
}

 4.特别的函数

在分子分母里同时乘(x-xi)

本条在之后埃尔米特插值法里会见到。

 5.优缺点

优点:结构紧凑、理论分析方便

缺点:改变一个节点,所有的差值基函数都会改变(基函数没有承袭性,于是引入牛顿差值)

(三)牛顿插值

0.前言

由于拉格朗日插值基函数不具有承袭性,于是思考是否能找到一组新的基函数,使得新增的基函数对原有函数进行修正,就能防止有新节点插入时Lagrange的重复计算,且前面的计算依然可以使用。

1.插值基函数的优化

首先作为插值基函数,要求“线性无关”。而任何一个n次多项式都可以表示为下图所示(n+1)个多项式的线性组合,并且他们线性无关,符合基函数要求。

 那么插值多项式:

 2.差商----各项系数

牛顿插值法的特点在于:每增加一个点,不会导致之前的重新计算,只需要算和新增点有关的就可以了。

公式推导如图所示:

 一二三...根据前一阶差商以此类推,我们可以画一个差商表:

 当我们新增一个点时,我们只需要在之前算得的基础上新增一项。

牛顿插值多项式的各项系数就是差分表中每列的第一项

3.代码的实现

//牛顿插值法,index表示插值点  
double Newton(int n, double x[], double y[], double index)
{
	double res = 0, sum;
	double temp = 1;//temp用来记录每一项的分母 
	int i, j;
	for (i = 0; i < n - 1; i++)
	{
		for (j = n - 1; j > i; j--)//计算差分 
		{
			temp = x[j] - x[j - i - 1];
			y[j] -= y[j - 1];
			y[j] /= temp;
		}
		sum = y[i + 1];//系数 

		for (int k = 0; k <= i; k++)//系数*插值基函数 (index代入公式中的x) 
		{
			sum *= (index - x[k]);
		}
		res += sum;
	}
	res += y[0];//y[0]就是我们的a[0] 
	return res;
}

4.优缺点

 优点:每增加一个点,不会导致之前的重新计算,只需要算和新增点有关的就可以了。

缺点:牛顿插值绘制的曲线与拉格朗日插值的曲线略有不同,但次数较高时,牛顿插值也会产生剧烈的震荡。另外,牛顿插值也存在连接处不光滑的缺陷。

(四) 埃尔米特插值

1.简介

埃尔米特插值在给定的节点处,不但要求 插值 多项式 的函数值与原函数值相同。 同时还要求在节点处,插值多项式的一阶直至指定阶的导数值,也与被插函数的相应阶导数值相等,这样的插值称为 埃尔米特 (Hermite)插值 。

n + 1个含函数值和导数值的插值条件可构造次数 ⩽ n 的 Hermite 插值多项式。

(下图来自https://blog.csdn.net/SanyHo/article/details/106849323

 2.基函数的确定

埃尔米特插值法有两个基函数,一个用于控制函数值,一个用于控制一阶导数值。

①控制函数值的基函数A(满足下图中紫色的三个条件)

 ②控制一阶导数值的基函数B(满足下图中紫色的三个条件)

 最后得出Hermite插值多项式H ( x )或者称为带倒数的插值多项式。其中l(小写的L) 就是拉格朗日插值的基函数。

3.代码

//埃尔米特插值
double Hermiter(int n, double x[3], double y[3], double y1[3], double index)
{
	double l[3], l1[3], a[3], b[3];
	for (int i = 0; i < 3; i++)
		l[i] = 1.0, l1[i] = 0.0;
	double res = 0;
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			if (j != i)
			{
				l1[i] += 1 / (x[i] - x[j]);
				l[i] *= (index - x[j]) / (x[i] - x[j]);
			}
		}
	}
	for (int i = 0; i < n; i++)
	{
		a[i] = (1 - 2 * (index - x[i]) * l1[i]) * (l[i] * l[i]);	//y[i]的系数a[i]
		b[i] = (index - x[i]) * (l[i] * l[i]);	//y[i]的导数的系数 b[i]
		res += (y[i] * a[i] + y1[i] * b[i]);
	}
	return res;
}

4.优缺点

埃尔米特插值得到的函数曲线更为平滑且接近原曲线,适用于除了拟合函数需要通过给定的点,还要求拟合的函数在给定点的导数也等于原函数在该点的导数。

不过也因此hermite依赖的更多,需要知道导数值。

(五)西楚霸王的自编之范德蒙德待定系数法

(本条是好哥们ZSD自己的思路,做了微微改动)

1.原理

原理来源于对线性代数的学习思考。运用到的知识主要是矩阵乘法和矩阵求逆。

文字难以表达,以三个点为例做如下图:

 2.代码

double Vandermonde(int n, double x[], double y[], double index)
{
	double a[3] = { 0 }, sum = 0;//a用于存放系数,sum用于累加y 
	double ni_vand[3][3], vand[3][3];//vand存范德蒙德,ni_vand存矩阵逆 
	int i, j;
	//初始化矩阵vand为范德蒙德行列式的形式 
	for (i = 0; i <= n; i++) 
	{
		for (j = 0; j <= n; j++) 
		{
			vand[i][j] = pow(x[j], i);
		}
	}
	double A = getA(vand, n + 1);//记录矩阵vand行列式的值
	getAStart(vand, ni_vand,n+1);//将vand矩阵的伴随矩阵存入ni_vand 

	//除以行列式得到vand的逆  
	for (i = 0; i <= n; i++) 
	{
		for (j = 0; j <= n; j++) 
        {
			ni_vand[i][j] /= A;
		}
	}//计算系数矩阵的逆 

	ChengJuChen(n,y,ni_vand,a);//矩阵乘法求系数a 

	sum += a[0];
	sum += a[1] * index;
	sum += a[2] * index * index;
	return sum;//y的多项式累加 
}

double getA(double a[3][3], int n)//计算|A|
{
	double A = 0;
	if (n == 2) {//二阶行列式 
		A += a[0][0] * a[1][1];
		A -= a[0][1] * a[1][0];
	}
	else if (n == 3) {//三阶行列式 
		A += a[0][0] * a[1][1] * a[2][2];
		A += a[0][1] * a[1][2] * a[2][0];
		A += a[0][2] * a[1][0] * a[2][1];
		A -= a[0][2] * a[1][1] * a[2][0];
		A -= a[0][1] * a[1][0] * a[2][2];
		A -= a[0][0] * a[1][2] * a[2][1];
	}
	return A;
}
void getAStart(double x1[3][3], double x2[3][3],int n)//计算A*(A的伴随矩阵)
{
	if (n == 2) //二阶矩阵直接运算
	{ 
		x2[0][0] = x1[1][1];
		x2[0][1] = -1 * x1[0][1];
		x2[1][0] = -1 * x1[1][0];
		x2[1][1] = x1[0][0];
		return;
	}
	if (n == 1)//一阶矩阵 
	{
		x2[0][0] = 1;
		return;
	}
	double temp[3][3];//中间变量数组 
	int i, j, k, t;
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			for (k = 0; k < n - 1; k++)
			{
				for (t = 0; t < n - 1; t++)
				{
					temp[k][t] = x1[k >= i ? k + 1 : k][t >= j ? t + 1 : t];//存放余子式矩阵 
				}
			}
			x2[j][i] = getA(temp, n - 1);//计算余子式的值 
			if ((i + j) % 2 == 1)//余子式的符号 
			{
				x2[j][i] = -x2[j][i];
			}
		}
	}
}


void ChengJuChen(int n, double x[3], double y[3][3], double c[3])//矩阵乘法
{
	int i, j;
	for (i = 0; i <=n; i++)
	{
		c[i] = 0;
		for (j = 0; j <= n; j++)
			c[i] += x[j]*y[j][i];
	}
}

3.优缺点

优点是很好理解,只用到了简单的矩阵运算,缺点和拉格朗日插值、牛顿插值类似,不再赘述。 

三、实验结果

1.主函数

实验一提供了三个点,所以调用两个点的一次插值法可以有三种不同的调用 。

int main()
{
	 double x[3] = { 1,2,3 };
	 double y[3] = { 0.367879441,0.135335283,0.049787068 };//老师给定三组数值
	 double y1[3] = { -0.367879441,-0.135335283,-0.049787068 };//导数
	 double x13[2] = { 1,3 };
	 double y13[2] = { 0.367879441,0.049787068 };//用于一次插值的(x1,y1)和(x3,y3) 
	 double y131[2] = { -0.367879441,-0.049787068 };
	//提供了三个点
	//调用两个点的一次插值法可以有三种不同的调用 

	//用前两个点插值
	printf("*********************************\n前两个点进行插值:\n");
	printf("拉格朗日一次插值:   %f\n", Lagrange(2, x, y, 2.6));
	printf("牛顿法一次插值:     %f\n", Newton(2, x, y, 2.6));
	//牛顿插值法会改变原数组,故调用完后重置 
	y[0] = 0.367879441, y[1] = 0.135335283, y[2] = 0.049787068;
	printf("埃尔米特插值法一次插值: %f\n", Hermiter(2, x, y,y1,2.6));
	printf("范德蒙德矩阵法一次插值: %f\n", Vandermonde(1,x,y,2.6));
	

	//用后两个点插值
	printf("*********************************\n后两个点进行插值:\n");
	printf("拉格朗日一次插值:   %f\n", Lagrange(2, x + 1, y + 1, 2.6));
	printf("牛顿法一次插值:     %f\n", Newton(2, x + 1, y + 1, 2.6));
	//重置原数组,理由同上 
	y[0] = 0.367879441, y[1] = 0.135335283, y[2] = 0.049787068;
	printf("埃尔米特插值法一次插值: %f\n", Hermiter(2, x+1, y+1, y1+1, 2.6));
	printf("范德蒙德矩阵法一次插值: %f\n", Vandermonde(1,x+1,y+1,2.6));
	
	//用第一个和第三个数据插值 
	printf("*********************************\n第一个和第三个数据插值:\n");
	printf("拉格朗日一次插值:   %f\n", Lagrange(2, x13, y13, 2.6));
	printf("牛顿法一次插值:     %f\n", Newton(2, x13, y13, 2.6));
	y[0] = 0.367879441, y[1] = 0.135335283, y[2] = 0.049787068;
	printf("埃尔米特插值法一次插值: %f\n", Hermiter(2, x, y13, y131, 2.6));
	printf("范德蒙德矩阵法一次插值: %f\n", Vandermonde(1,x13,y13,2.6));
	

	//用三个数进行插值
	printf("*********************************\n三个点进行插值:\n");
	printf("拉格朗日二次插值:   %f\n", Lagrange(3, x, y, 2.6));
	printf("牛顿法二次插值:     %f\n", Newton(3, x, y, 2.6));
	y[0] = 0.367879441, y[1] = 0.135335283, y[2] = 0.049787068;
    y13[0] = 0.367879441,y13[1]=0.049787068 ;
	printf("埃尔米特插值法二次插值: %f\n", Hermiter(3, x, y, y1, 2.6));
	printf("范德蒙德矩阵法二次插值: %f\n", Vandermonde(2,x,y,2.6));
	return 0;
}

2.运行结果

可以看到埃尔米特的精确度确实更高。而其它三种方法本质是差不多的,所以结果也呈现得相似。

3.四种方法的比较

①拉格朗日插值法:结构紧凑、理论分析方便改变一个节点,所有的差值基函数都会改变(基函数没有承袭性)。时间复杂度:O(n²)

②牛顿插值法:每增加一个点,不会导致之前的重新计算,只需要算和新增点有关的就可以了。但是连接处不光滑数较高时,牛顿插值也会产生剧烈的震荡。时间复杂度:O(n²)

③埃尔米特插值法:准确度更高,拟合曲线更平滑,但依赖更多。时间复杂度:O(n²)

④范德蒙德待定系数法:方便理解,但计算复杂。移植性最强。时间复杂度:O(n⁴)

四、实验感受

这次实验让我感受到了数学的压迫QAQ,查阅很多资料去学习原理和编程相比是一个更漫长的过程,花费很多时间,但同时也有很多收获。在证明方面还有一些囫囵吞枣的bug,希望有所改进。

  • 5
    点赞
  • 17
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

程序和三三总有一个能跑

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值