[数值分析]插值 - 牛顿插值法

clc;clear all;close all;

n = 1:1:10000;
sig = sin(n/100);

k = 0;
for i = 1:4:10000
    k = k + 1;
    x_arr = n(i:i+3);
    y_arr = sig(i:i+3);
    arr_new_y(k) = NewtonInterpolation(x_arr,y_arr,4,mean(x_arr));
    arr_new_x(k) = mean(x_arr);
end

figure(1);
plot(n,sig);hold on;
scatter(arr_new_x,arr_new_y);grid on;

function [y_return] = NewtonInterpolation(x_arr,y_arr,n,x_in)

    % 提早返回
    if(n <= 1) 
        y_return = 0;
        return;
    end

    % 商差缓存申请与计算
    buffer_len   = (n * (n+1)) / 2;
    quo_diff_arr = zeros(buffer_len,1);
    for i = 1:1:n
        quo_diff_arr(i) = y_arr(i); % 三角第一列为原始y值用于迭代
    end
    
    % 开始计算插值
    % x y
    % * *
    % * **
    % * ***
    % * ****
    
    y_arr_index     = 0; % 三角前一排y值
    y_new_arr_index = n; % 三角后一排y值 提前赋值
    
    for diff_len = 1:1:n-1 % 差商距离 如果n为4 则距离由1-3
        for index = 1:1:n-diff_len % x值按差商距离扫描一次商差计算
            y_arr_index     = y_arr_index + 1;
            y_new_arr_index = y_new_arr_index + 1;

            % 计算差商
            quo_diff_arr(y_new_arr_index) = ...
                (quo_diff_arr(y_arr_index+1) - quo_diff_arr(y_arr_index))/...
                (x_arr(index+diff_len) - x_arr(index));
            
        end
        y_arr_index = y_arr_index + 1;
    end
  
    % 代入多项式 f(x) = f(x0) + f(x1)(x-x0) + f(x2)(x-x1)(x-x0) + ...+ f(xn)(x-xn-1)~(x-x0);
    sum     = 0;
    step    = n;
    index   = 1;
    for i = 1:1:n
        j   = i - 1;
        mul = 1;
        while(j > 0)
            mul = mul*(x_in - x_arr(j));
            j = j - 1;
        end
        sum = sum + quo_diff_arr(index)*mul;
        
        index   = index + step;
        step    = step - 1;    
    end
    
    % 返回结果
    y_return = sum;
    return;
end