文章内容:使用三次多项式对y = sin2πx进行插值(确切说应该叫拟合,下面姑且还叫插值,因为题目就是这样出的)。下面只涉及解题思路,不涉及相关概念。
一、解题思路
二、数据处理
三、运行结果
1.三次多项式系数
2.插值结果与y = sin2πx的比较
四、代码实现(matlab)
%创作者:壹大子皿
%时间:2019/12/18
%地点:北京
%****************************************
%main.m
%****************************************
clear,clc
x0 = 0:0.05:1;
y0 = sin(2*pi*x0);
for i = 1:length(y0)
if abs(y0(i)) < 1e-15
y0(i) = 0;
end
end
w = polyReg(x0', y0', 3);
w0 = w(4), w1 = w(1), w2 = w(2), w3 = w(3)
x = 0:0.01:1;
x = x';
y = [x x.^2 x.^3 ones(length(x),1)] * w;
%绘制结果
plot(x, sin(2*pi*x), 'b-')
hold on
plot(x, y, 'r-.')
legend('y = sin2πx','y = w0 + w1*x + w2*x^2 + w3*x^3')
xlabel('x轴'),ylabel('y轴')
title('对 y = sin2πx 进行3次多项式插值')
%****************************************
%polyReg.m
%****************************************
%多项式线性回归
function W = polyReg(x0, y0, n)
%线性回归,最小二乘法
%(x0,y0):已有数据集,均为m*1
%n:多项式阶数
m = length(x0); %数据集大小
X = zeros(m, n+1);
for i = 1:n
X(:,i) = x0 .^ i;
end
X(:,n+1) = 1;
T = X'* X;
[~,s]=eig(T);
flag = true;
for i = 1:length(s)
if s(i,i) <= 0
flag = false;
end
end
if rank(T) == n+1 || flag == true %如果T为满秩矩阵或正定矩阵
W = T \ X' * y0;
end
end