【数值分析】拉格朗日插值与牛顿插值

在工程应用和科学研究中，经常要研究变量之间的关系y=f(x)。但对于函数f(x)，常常得不到一个具体的解析表达式，它可能是通过观测或实验得到的一组数据(x,f(x))，x为一向量;或则是解析表达式非常复杂，不便于计算和使用。因此我们需要寻找一个计算比较简单的函数S(x)近似代替f(x)，并使得S(x)=f(x)，这种方法就称为插值法。

常用的插值法有：
一维插值法：拉格朗日插值、牛顿插值、分段低次插值、埃尔米特插值、样条插值。
二维插值法：双线性插值、双二次插值。

拉格朗日插值法

已知函数f(x)的n+1个互异的节点

处的函数值

，则其拉格朗日插值多项式可以写为：



其中，

为插值基函数，其表达式为：

牛顿插值法

已知函数f(x)的n+1个互异的节点

处的函数值

，则其牛顿插值多项式可以写为：



其中，

为f(x)的k阶差商（也叫均差），可以表示如下：



也可以由函数值

线性表示为：



根据上述基本原理和公式，很容易编程实现。我们假设根据下面的数据表，来分别用拉格朗日插值和牛顿插值来计算f(8.4)的近似值：

x	8.1	8.3	8.6	8.7
f(x)	16.94410	17.56492	18.50515	18.82091

具体代码实现如下：

#include<iostream>
#include<vector>
using namespace std;
//-----------------拉格朗日插值法BEGIN---------------------//
double Lagrange(vector<double> x,vector<double> y ,double X)//x,y分别为x和f(x)的值，X为要求的点，返回值为f(X)
{
double result=0;
double temp=1;
for(int i=0;i<x.size();i++)
{
temp=1;
for(int j=0;j<x.size();j++)
{
if(j!=i)
{
temp=temp*(X-x.at(j))/(x.at(i)-x.at(j));
}
}
result+=temp*y.at(i);
}
return result;
}
//-----------------拉格朗日插值法END---------------------//
//-----------------牛顿法BEGIN---------------------//
double DifferenceQuotient(vector<double> x,vector<double> y ,int k)//计算差商
{
double result=0;
for(int i=0;i<=k;i++)
{
double temp=1;
for(int j=0;j<=k;j++)
{
if(i!=j)
{
temp=temp/(x.at(i)-x.at(j));
}
}
temp=y.at(i)*temp;
result+=temp;
}
return result;
}
double Newton(vector<double> x,vector<double> y ,double X)
{
double result=y.at(0);
double temp=1;
for(int i=1;i<x.size();i++)
{
temp=1;
for(int j=0;j<i;j++)
{
temp*=(X-x.at(j));
}
result+=DifferenceQuotient(x,y,i)*temp;
}
return result;
}
//-----------------牛顿法END---------------------//
void main()
{
vector<double> x;
vector<double> y;
//这里输入x的值,这里使用向量vector是为了方便添加数据点，可以根据实际的观测点更改
x.push_back(8.1);
x.push_back(8.3);
x.push_back(8.6);
x.push_back(8.7);

//这里输入f(x)的值
y.push_back(16.94410);
y.push_back(17.56492);
y.push_back(18.50515);
y.push_back(18.82091);

cout.precision(10);//设置显示精度
//下面是根据上面的4个样本点及其函数值来分别使用两种插值法计算在x=8.4处的函数值
cout<<"使用拉格朗日插值法：";
cout<<Lagrange(x,y,8.4)<<endl;
cout<<"使用牛顿插值法：";
cout<<Newton(x,y,8.4)<<endl;
}

程序运行结果如下：



优缺点比较：
拉格朗日插值法：插值多项式和插值基函数的形式对称，容易编程。但是，增加节点时，需要重新计算每一个插值基函数。
牛顿插值法：当插值节点增加时，之前已计算的结果仍然能用，每增加一个节点，只要再增加一项即可，从而避免了重复性计算。

Matlab实现多种插值函数

现在也有很多人使用Matlab来进行算法的仿真，我在这里把大二时数学建模整理的插值算法函数也共享出来，链接为：http://download.csdn.net/detail/tengweitw/8387451具体的使用说明在文件中都有说明。下面就拿我们刚才所讲的拉格朗日插值和牛顿插值来举例说明，还是使用上面的数据表，则拉格朗日插值函数如下：

function f = Language(x,y,x0)
%x y为坐标向量  x0为插值点的x坐标|| f0为x0对应的值

syms t;
if(length(x) == length(y))
n = length(x);
else
disp('x和y的维数不相等！');
return;
end                                      %检错

f = 0.0;
for(i = 1:n)
l = y(i);
for(j = 1:i-1)
l = l*(t-x(j))/(x(i)-x(j));
end;
for(j = i+1:n)
l = l*(t-x(j))/(x(i)-x(j));      %计算拉格朗日基函数
end;

f = f + l;                           %计算拉格朗日插值函数
simplify(f);                         %化简

if(i==n)
if(nargin == 3)
f = subs(f,'t',x0);          %计算插值点的函数值
else
f = collect(f);          %将插值多项式展开
f = vpa(f,6);                %将插值多项式的系数化成6位精度的小数
end
end
end

牛顿插值函数如下：

function f = Newton(x,y,x0)
%x y为坐标向量  x0为插值点的x坐标|| f0为x0对应的值
syms t;

if(length(x) == length(y))
n = length(x);
c(1:n) = 0.0;
else
disp('x和y的维数不相等！');
return;
end

f = y(1);
y1 = 0;
l  = 1;

for(i=1:n-1)
for(j=i+1:n)
y1(j) = (y(j)-y(i))/(x(j)-x(i));
end
c(i) = y1(i+1);
l = l*(t-x(i));
f = f + c(i)*l;
simplify(f);
y = y1;

if(i==n-1)
if(nargin == 3)
f = subs(f,'t',x0);
else
f = collect(f);                %将插值多项式展开
f = vpa(f, 6);
end
end
end

为了使用上面两个函数，脚本文件如下：

clear all
clc
format long
format compact

x=[8.1 8.3 8.6 8.7 ];
y=[ 16.94410 17.56492 18.50515 18.82091];
x0=8.4;
disp('拉格朗日插值法：')
disp(Language(x,y,x0))
disp('牛顿插值法：')
disp(Newton(x,y,x0))

结果显示如下：



原文：http://blog.csdn.net/tengweitw/article/details/43025225

作者：nineheadedbird