牛頓插值 C++ 和 Matlab實現

牛頓插值法

概述

均差(差商)

先寫一下均差定義
點 x0和xk

x_0 和 x_k

的一階均差 f[x0,xk]=f(xk)−f(x0)xk−x0

f[x_0,x_k] = \frac{f(x_k) - f(x_0)}{x_k- x_0}

利用上式 關於 x0,x1,xk

x_0 , x_1 , x_k

二階均差 為 f[x0,x1,xk]=f[x0,xk]−f[x0,x1]xk−x1

f[x_0,x_1,x_k] = \frac{f[x_0,x_k] - f[x_0,x_1]}{x_k- x_1}

依次類推得 k 階均差

f[x0,x1,...,xk]=f[x0,...,xk−2,xk]−f[x0,x1,...,xk−1]xk−xxk−1

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}}

均差(差商)表

注意每列的左右關係

牛頓插值

給出兩點及函式值 x0,f(x0),x1,f(x1)

x_0,f(x_0),x_1,f(x_1)

插值多項式代表直線點斜式

f(x1)=f(x0)+f(x1)−f(x0)x1−x0×(x1−x0)

f(x_1) = f(x_0) + \frac{f(x_1) - f(x_0)}{x_1 - x_0} \times (x_1 - x_0)

用均差替換

f(x1)=f(x0)+f[x0,x1]×(x1−x0)

f(x_1) = f(x_0) + f[x_0,x_1] \times (x_1 - x_0)

三個節點 x0,f(x0),x1,f(x1),x2,f(x2)

x_0,f(x_0),x_1,f(x_1),x_2,f(x_2)

插值多項式

f(x2)=f(x0)+f[x0,x1]×(x2−x0)+f[x0,x1,x2]×(x2−x1)×(x2−x0)

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 個點x0,f(x0),x1,f(x1),...,xn−1,f(xn−1),x,f(x)

x_0,f(x_0),x_1,f(x_1),...,x_{n-1},f(x_{n-1}),x,f(x)

f(x)=f(x0)+f[x0,x1]×(x−x0)+f[x0,x1,...,xn−1,x]×(x−xn−1)×...×(x−x1)×(x−x0)

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
                    Copyright 2017 Gy 
tips : Gy是高玉的標誌，可見logo
/////////////////////////////////////////////////////////////////
*/ 

#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)