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
[数值分析]插值 - 待定系数法
最新推荐文章于 2023-02-04 10:37:04 发布