在数学建模中,我们拿到的数据经常会有缺失值,在缺失值不是很多的情况下,我们在数据预处理阶段会采用插值方法来将数据补齐,之后再开始我们的建模。
目录
1.Matlab 实现分段线性插值
2.拉格朗日插值多项式
3.牛顿(Newton)插值
编辑 4.埃尔米特(Hermite)插值
1.Matlab 实现分段线性插值
Matlab 中有现成的一维插值函数 interp1。
y=interp1(x0,y0,x,'method')
method 指定插值的方法,默认为线性插值。其值可为: 'nearest' 最近项插值, 'linear' 线性插值 'spline' 逐段 3 次样条插值, 'cubic' 保凹凸性 3 次插值。
所有的插值方法要求 x0 是单调的。
示例:
% 已知的点(x坐标和对应的y坐标) x = [1, 3, 5, 7, 9]; y = [2, 4, 6, 8, 10]; xi = 1:0.5:9; % 从1到9,步长为0.5 % 使用interp1函数进行线性插值 yi = interp1(x, y, xi, 'linear'); plot(x, y, 'o', xi, yi, '-'); legend('原始点', '线性插值'); xlabel('x'); ylabel('y'); title('线性插值示例'); grid on;disp('原始点y:'); disp(y); disp('插值后yi:'); disp(yi);
2.拉格朗日插值多项式
首先构造一组基函数:
上式称为n 次 Lagrange 插值多项式。
Matlab 中没有现成的 Lagrange 插值函数,必须编写一个 M 文件实现 Lagrange 插值。 设n 个节点数据以数组 x0, y0输入(注意 Matlat 的数组下标从 1 开始),m 个插值 点以数组 x 输入,输出数组 y 为m 个插值。编写一个名为 lagrange.m 的 M 文件:
function y=lagrange(x0,y0,x); n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end
代入数据:
x0 = [1, 3, 5, 7, 9]; y0 = [2, 4, 6, 8, 10]; x = 1:0.5:9;y=lagrange(x0,y0,x); disp('原始点y0:'); disp(y0); disp('插值后y:'); disp(y);
输出结果:
3.牛顿(Newton)插值
运行程序:
function yi = newtonInterpolation(x, y, xi) n = length(x); % 已知数据点的数量 if n ~= length(y) error('x和y的长度必须相等'); end % 构建差分表 diffTable = zeros(n, n); diffTable(:, 1) = y(:); % 第一列是y值 for j = 2:n for i = 1:n-j+1 diffTable(i, j) = (diffTable(i+1, j-1) - diffTable(i, j-1)) / (x(i+j-1) - x(i)); end end if ~isvector(xi) yi = diffTable(1, 1); % 插值多项式的常数项 for k = 2:n p = 1; for j = 1:k-1 p = p * (xi - x(j)); end yi = yi + diffTable(1, k) * p; end else yi = zeros(size(xi)); for i = 1:length(xi) yi_temp = diffTable(1, 1); % 插值多项式的常数项 for k = 2:n p = 1; for j = 1:k-1 p = p * (xi(i) - x(j)); end yi_temp = yi_temp + diffTable(1, k) * p; end yi(i) = yi_temp; end end end
代入数据:
x = [0, 1, 2, 3, 4]; y = [1, 0.8, 0.9, 0.1, -0.8]; xi =0:0.5:5; yi = newtonInterpolation(x, y, xi); % 绘图 plot(x, y, 'o', xi, yi, '-'); legend('原始数据点', '牛顿插值'); xlabel('x'); ylabel('y'); title('牛顿插值示例'); grid on;disp('原始点y:'); disp(y); disp('插值后yi:'); disp(yi);
运行结果:
4.埃尔米特(Hermite)插值
如果对插值函数,不仅要求它在节点处与函数同值,而且要求它与函数有相同的一 阶、二阶甚至更高阶的导数值,这就是 Hermite 插值问题。
运行程序:
function y=hermite(x0,y0,y1,x); n=length(x0);m=length(x); for k=1:m yy=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j~=i h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2; a=1/(x0(i)-x0(j))+a; end end yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i)); end y(k)=yy; end
代入数据:
x0=[0.1, 0.2,0.3,0.4,0.5]; y0=[1,2,3,4,5]; y1=[2,4,6,8,10];xi =0:0.5:5; yi = hermite(x0, y0,y1, xi); % 绘图 plot(x0, y0, 'o', xi, yi, '-'); legend('原始数据点', 'hermite插值'); xlabel('x'); ylabel('y'); grid on;disp('原始点y0:'); disp(y0); disp('插值后yi:'); disp(yi);
运行结果: