[数值分析]插值 - 待定系数法

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) = MUCInterpolation(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] = MUCInterpolation(x_arr,y_arr,n,x_in)

    % 提早返回
    if(n <= 1) 
        y_return = 0;
        return;
    end
    
    % 计算公式 解矩阵 ax0^(n-1) + bx0^(n-2) + ... + kx0 + c = y(x0);
    % x0 遍历 x1 ~ xn输入
    
    % 计算x轴相对偏移 防止矩阵解计算量过大溢出
    pos_fix = x_arr(1);
    for i = 1:1:n
        x_arr(i) = x_arr(i) - pos_fix;
    end
    x_in = x_in - pos_fix;
    
    % 构建矩阵 解待定系数
    mat = zeros(n,n);
  
    for i = 1:1:n % 行循环
        for j = 1:1:n % 列循环
            mul = 1;
            % 常数项
            if(j == n)
                mat(i,j) = 1;
            else
                for k = 1:1:n - j % 系数循环
                    mul = mul*(x_arr(i));
                end
                mat(i,j) = mul;
            end
        end
    end
    
    % 解矩阵
    solve_arr = mat\y_arr';
    
    % 通过系数计算新y值
    sum = 0;
    for i = 1:1:n % 行循环
        mul = 1;
        if(i == n)
            sum = sum + solve_arr(n);
        else
            for j = 1:1:n - i % 行循环
                mul = mul * x_in;
            end
            sum = sum + solve_arr(i)*mul;
        end
    end
    
    % 返回结果
    y_return = sum;
    return;
end