目录
一、实验内容
(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,希望有所改进。