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

#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实现多种插值函数

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