牛顿插值 C++ 和 Matlab实现

牛顿插值 C++ 和 Matlab实现

牛顿插值法

概述

均差(差商)

先写一下均差定义
点 $x_0 和 x_k $ 的一阶均差 $ f[x_0,x_k] = \frac{f(x_k) - f(x_0)}{x_k- x_0}$
利用上式 关于 $x_0 , x_1 , x_k $二阶均差 为 $ f[x_0,x_1,x_k] = \frac{f[x_0,x_k] - f[x_0,x_1]}{x_k- x_1} $
依次类推得 k 阶均差
$$f[x_0,x_1,...,x_k]=\frac{f[x_0,...,x_{k-2},x_k]-f[x_0,x_1,...,x_{k-1}]}{x_k-x_{x_k-1}}$$
均差(差商)表

注意每列的左右关系

牛顿插值

给出两点及函数值 $x_0,f(x_0),x_1,f(x_1)$
插值多项式代表直线点斜式 $$f(x_1)=f(x_0)+\frac{f(x_1)-f(x_0)}{x_1-x_0}\times (x_1-x_0)$$
用均差替换
$$f(x_1)=f(x_0)+f[x_0,x_1]\times (x_1-x_0)$$
三个节点 $x_0,f(x_0),x_1,f(x_1),x_2,f(x_2)$
插值多项式
$$f(x_2)=f(x_0)+f[x_0,x_1]\times (x_2-x_0)+f[x_0,x_1,x_2]\times (x_2-x_1)\times (x_2-x_0)$$
依次类推
n 个点$x_0,f(x_0),x_1,f(x_1),...,x_{n-1},f(x_{n-1}),x,f(x)$

$$f(x)=f(x_0)+f[x_0,x_1]\times (x-x_0)+f[x_0,x_1,...,x_{n-1},x]\times (x-x_{n-1})\times ...\times (x-x_1)\times (x-x_0)$$

说明

上述公式是我用mathjax写的，如有错误请联系我修正
敬请指正
概述省略了大量推导过程，请查阅详细推导资料

代码

C++ 实现

结果和输入格式在注释中

/*
/
                              结果 
/
输入需要插值的数目 : 6
Xi : 0.4 0.55 0.65 0.8 0.9 1.05
Yi : 0.41075 0.57815 0.69675 0.88811 1.02652 1.25382
牛顿插值均差表 :
   0.41075          0          0          0          0          0
   0.57815      1.116          0          0          0          0
   0.69675      1.186       0.28          0          0          0
   0.88811    1.27573   0.358933   0.197333          0          0
   1.02652     1.3841   0.433467   0.212952  0.0312381          0
   1.25382    1.51533   0.524933   0.228667  0.0314286 0.00029304
输入需要插值的x :
0.596
y : 0.631917
/
tips : 运行环境 DevC++ 5.9.2
/
*/ 

#include<iostream> 
#include<iomanip>
#include<cstring>
#define maxn 100
using namespace std;

double x[maxn];
double matrix[maxn][maxn];
//存放均差表矩阵 
int n;

double newton(double & y1,double & y2,double & x1,double & x2);
//插值计算表达式 
int mat_disp();
//矩阵现实 
 
int main()
{
	double ans,xi,t;  
	cout << "输入需要插值的数目 : " ; 
	cin >> n;
	memset(matrix,0,sizeof(matrix));
	cout << "Xi : ";
	for(int i = 0;i < n;i++)
		cin >> x[i];
	cout << "Yi : ";
	for(int i = 0;i < n;i++)
		cin >> matrix[i][0];
	for(int i = 1;i < n;i++)
		for(int j = i;j < n;j++)
		{
			matrix[j][i] = newton(matrix[j][i-1],matrix[j-1][i-1],x[j],x[j-i]);
		}
	mat_disp();
	cout << "输入需要插值的x : " << endl;
	cin >> xi;
	ans = matrix[0][0];
	for(int i = 1;i < n;i++)
	{
		t = matrix[i][i];
		for(int j = i-1;j >= 0;j--)
			t *= (xi - x[j]);
		ans += t;
	}
	cout << "y : " << ans << endl;
	return 0;
}

int mat_disp()
{
	cout << "牛顿插值均差表 : " << endl;
	for(int i = 0;i < n;i++)
	{
		for(int j = 0;j < n;j++)
		{
			cout << setw(10) << matrix[i][j] << ' ';
		}
		cout << endl;
	}
	return 0;
}

double newton(double & y1,double & y2,double & x1,double & x2)
{
	return (y2-y1)/(x2-x1); 
}

Matlab 实现

将该函数存为m文件

function res = Newton(x,y,xi)
    xl = length(x);
    yl = length(y);
    if(xl ~= yl)
        disp('向量长度不等')
        return
    end
    m = zeros(xl);
    for i = 1:xl
        m(i,1) = y(i);
    end
    for i = 2:xl
        for j = i:xl
            m(j,i) = ((m(j,i-1) - m(j-1,i-1))/(x(j) - x(j-i+1)));
        end
    end
    %disp(m);
    res = m(1,1);
    for i = 2:xl
        t = m(i,i);
        for j = i-1:-1:1
            t = t * (xi - x(j));
        end
        disp(res);
        res = (res + t);
        disp(res);
    end
end

调用下面语句测试函数

x = [0.4,0.55,0.65,0.8,0.9,1.05]
y = [0.41075,0.57815,0.69675,0.88811,1.02652,1.25382]
xi = 0.596
Newton(x,y,xi)