拉格朗日插值

数学原理

此处的拉格朗日插值均为多项式插值,固定下节点,多项式立刻确定下来。

n+1个互异节点满足插值条件的n次拉格朗日插值多项式为:

插值 的权重 如果处理_插值


插值 的权重 如果处理_i++_02


插值 的权重 如果处理_ci_03

代码实现

算法实现过程:
1.获取节点个数
2.将节点的x和y存入两个数组中
3.在进行拉格朗日插值中,利用两层循环,外循环累加每一项lk(x)的的值,内循环计算各项lk(x)的值(i!=j时进行累乘)

#include<iostream>
#include<stdlib.h>
using namespace std;

float x[10];
float y[10];
int count;

void load_data()
{
	cout<<"the count is";
	cin>>count;
	for(int i=0;i<count;i++)
	{
		cin>>x[i]>>y[i];
		cout<<"\n";
	}
}

void print_data()
{
	for(int i=0;i<count;i++)
	{
		cout<<"x is "<<x[i]<<" y is "<<y[i]<<'\n';
	}
}


float Lagrange_interpolation(float value_x)
{
	float sum_all=0;
	float sum_single;
	for (int i=0;i<count;i++)
	{
		sum_single=y[i];
		for(int j=0;j<count;j++)
		{
			if(i==j)
				continue;
			else
			{
				sum_single*=(value_x-x[j])/(x[i]-x[j]);
			}
		}
		//cout<<"the sum_single is "<<sum_single<<endl;
		sum_all+=sum_single;
	}
	return sum_all;
}

void main()
{
	float value_x;
	float result;
	load_data();
	print_data();
	cout<<"please input the value_x ";
	cin>>value_x;
	result=Lagrange_interpolation(value_x);
	cout<<"the output is "<<result;
	system("pause");
}

埃特金插值

数学原理

采用逐次线性插值的方法,可以灵活地增加插值节点,且具有所谓的“承袭性”

插值 的权重 如果处理_插值_04


上图即为原理图,每一次新增的插值都是采用线性插值,即取两个节点。

举例:

插值 的权重 如果处理_i++_05


得到的表达式将与拉格朗日插值相同。

特点:

1.将一个高次插值过程归结为线性插值的多次重复

2.埃特金插值表中的每个数据均可视为差值结果;这些数据的一致程度即可判断差值结果的精度

3.可以逐行生成插值表,每做一步便检查一下计算结果的精度,如不满足精度,则增加一个节点再算,直到满足精度为止

4.在插值节点较多的情况下,可以降低插值次数

代码实现
1.建立了数组node用于存放节点的x坐标,建立数组result用于存放不同阶数的插值,为一个三角矩阵
2.数组result的存放格式为:行代表的是插值的阶数,列代表的是对应阶数的不同的插值项
3.具体的计算逻辑见Aitken_interpolation(float x),主要是找前一阶的插值(首项以及对应项),前一阶的插值找到后就可以找到对应的节点值

#include <iostream>
#include <stdlib.h>
using namespace std;

int count;
float result[10][10];
float node[10];

void Load_data()
{
	cout<<"please input the count ";
	cin>>count;

	cout<<"input the node and result"<<endl;
	for(int i=0;i<count;i++)
	{
		cin>>node[i]>>result[0][i];
		cout<<endl;
	}
}

void Print_result()
{
	for (int i=0;i<count;i++)
	{
		for (int j=0;j<count;j++)
		{
			cout<<result[i][j]<<" ";
		}
		cout<<endl;
	}
}

void Aitken_interpolation(float x)
{
	for (int i=1;i<count;i++)
	{
		for (int j=i;j<count;j++)
		{
			result[i][j]=result[i-1][j]*(x-node[i-1])/(node[j]-node[i-1])+result[i-1][i-1]*(x-node[j])/(node[i-1]-node[j]);
		}
	}
}

void main()
{
	float x;
	Load_data();
	cout<<"please input the value_x ";
	cin>>x;
	Aitken_interpolation(x);
	Print_result();
	system("pause");
}

PS:如何实现尾部插入不需全部重新计算?

牛顿插值

数学原理

差商的定义:

一阶:

插值 的权重 如果处理_ci_06


二阶:

插值 的权重 如果处理_i++_07


n阶差商:

插值 的权重 如果处理_插值_08


差商表:

插值 的权重 如果处理_插值_09


建立差商表后,利用差商的定义式进行递推:

插值 的权重 如果处理_插值_10


插值 的权重 如果处理_插值_11


最后留出余项:

插值 的权重 如果处理_插值_12


可以发现,实际用到的各阶差商值为差商表中的对角线项。

插值 的权重 如果处理_插值_13


代码实现

1.建立两个数组,数组node代表节点的x值,二维数组quotient_table代表的是各阶的差商的表

2.差商表的每一行代表的是阶数,每一列代表的是同一阶数的不同项的差商

3.在计算差商时,需要找到上一阶的差商值(即上一阶的同列以及上一列),也要找到节点(一个节点的索引为该列序号,另一个节点索引应为该列序号减去阶数),具体见Build_quotient_table()

4.实现牛顿插值时,使用temp_x作为中间变量对(x-xi)的累乘进行保存,减少计算次数

#include <iostream>
#include <stdlib.h>
using namespace std;

int count;
float quotient_table[10][10];
float node[10];

void Load_Data()
{
	cout<<"please input the count ";
	cin>>count;
	cout<<"\nplease input the node and zero order of quotient"<<endl;
	for (int i=0;i<count;i++)
	{
		cin>>node[i]>>quotient_table[0][i];
		cout<<endl;
	}
}

void Print_quotient_table()
{
	for (int i=0;i<count;i++)
	{
		for (int j=0;j<count;j++)
		{
			cout<<quotient_table[i][j]<<" ";
		}
		cout<<endl;
	}
}

/*
 *@name Build_quotient_table:建立差商表
 *@return none
*/
void Build_quotient_table()
{
	for (int i=1;i<count;i++)
	{
		for (int j=i;j<count;j++)
		{
			//二维数组每一行i对应的是i阶差商,每一列j对应的是在该阶差商中的不同项之间的差商
			quotient_table[i][j]=(quotient_table[i-1][j]-quotient_table[i-1][j-1])/(node[j]-node[j-i]);
		}
	}
}

/*
 *@name Newton_interpolation:利用牛顿插值进行计算
 *@param x:变量x
 *@return:牛顿插值得到的结果
*/
float Newton_interpolation(float x)
{
	float sum=0;
	float temp_x=1;
	for (int i=0;i<count;i++)
	{
		sum+=temp_x*quotient_table[i][i];
		temp_x*=(x-node[i]);
	}
	return sum;
}

void main()
{
	Load_Data();
	Build_quotient_table();
	Print_quotient_table();
	cout<<"the result is "<<Newton_interpolation(0.596)<<endl;
	system("pause");
}