目录
线性插值
原理
流程图
单个点的线性插值代码
X=[0.2 0.4];
Y=[21 25];
x=0.7;
x0=X(1)
y0=Y(1);
x1=X(2);
y1=Y(2);
L0=(x-x1)/(x0-x1);
L1=(x-x0)/(x1-x0);
y=y0*L0+y1*L1;
多个点的线性插值代码
time = [1 3 8 12 15 20 24];
tem = [8 9 16 23 22 18 10];
time_i = 1:0.01:24;
tem_i = self_interp1(time,tem,time_i,'linear');
plot(time,tem,'o',time_i,tem_i);
%自己写一个interp1类似功能的接口
%在这个接口中参数x需要从大到小排列,y随意
function [yi]=self_interp1(x,y,xi,method)
% 初始化yi,给它xi对应的列
col_xi = size(xi,2);
yi = zeros(1,col_xi);
% 检测使用的插值方法 这里期望的是'linear'
if strcmp(method,'linear')
% 找到每个xi在x序列中的位置
col_x = size(x,2);
for i = 1:col_xi,
for j = 1:col_x-1,
% 假如需要计算插值公式
if x(j+1) > xi(i),
yi(i) = y(j)+(y(j+1)-y(j))/(x(j+1)-x(j))*(xi(i)-x(j));
break;
end
% 假如插值处的数据已经测得了,就直接把值给它,节约计算资源
if x(j) == xi(i),
yi(i) = y(j);
break;
end
end
% 以上没有把最后一个数据点考虑进去,需要加上
yi(col_xi) = y(col_x);
end
else
error('插值方法请选择(linear)\n');
end
end
抛物插值
原理
流程图
代码
X=[0.2 0.4 0.6];
Y=[21 25 23];
x0=X(1);
x1=X(2);
x2=X(3);
y0=Y(1);
y1=Y(2);
y2=Y(3);
x=0.7;
L0=(x-x1)*(x-x2)/(x0-x1)/(x0-x2);
L1=(x-x0)*(x-x2)/(x1-x0)/(x1-x2);
L2=(x-x0)*(x-x1)/(x2-x0)/(x2-x1);
y=y0*L0+y1*L1+y2*L2
拉格朗日插值
代码
x=[0.2 0.4 0.6 0.8];
y=[21 25 23 20];
yh=lagrange(x,y,0.7)
function yh=lagrange (x,y,xh)
n = length(x);
m = length(xh);
yh = zeros(1,m);
c1 = ones(n-1,1);
c2 = ones(1,m);
for i=1:n
xp = x([1:i-1 i+1:n]);
yh = yh + y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));
end
end
牛顿插值
原理
代码
xi=[1 4 9];
yi=[1 2 3];
x=7;
p= Newton_fun(x,xi,yi)
function p= Newton_fun(x,xi,yi)
n=length(xi);
f=zeros(n,n);
% 对差商表第一列赋值
for k=1:n
f(k)=yi(k);
end
% 求差商表
for i=2:n % 差商表从0阶开始;但是矩阵是从1维开始存储!!!!!!
for k=i:n
f(k,i)=(f(k,i-1)-f(k-1,i-1))/(xi(k)-xi(k+1-i));
end
end
disp('差商表如下:');
disp(f);
%求插值多项式
p=0;
for k=2:n
t=1;
for j=1:k-1
t=t*(x-xi(j));
disp(t)
end
p=f(k,k)*t+p;
disp(p)
end
p=f(1,1)+p;
end
分段线性插值
原理
代码
x = [1 3 8 12 15 20 24];
y = [8 9 16 23 22 18 10];
yy=fdxx(x,y,7)
function yy=fdxx(x,y,xx)
n=size(x,2);
for i=1:n-1
if x(i)<xx&&xx<x(i+1)
L1=(xx-x(i+1))/(x(i)-x(i+1));
L2=(xx-x(i))/(x(i+1)-x(i));
yy=L1*y(i)+L2*y(i+1);
break;
elseif x(i)==xx
yy=y(i);
end
end
end
分段抛物插值
原理
代码
x = [1 3 8 12 15 20 24];
y = [8 9 16 23 22 18 10];
y=fenduanpaowu(x,y,7)
function y=fenduanpaowu(xi,yi,x)
n=size(xi,2);
if x<xi(2)
L1=(x-xi(1))*(x-xi(3))/(xi(1)-xi(2))/(xi(1)-xi(3));
L2=(x-xi(1))*(x-xi(3))/(xi(2)-xi(1))/(xi(2)-xi(3));
L3=(x-xi(1))*(x-xi(2))/(xi(3)-xi(1))/(xi(3)-xi(2));
y=L1*yi(1)+L2*yi(2)+L3*yi(3);
elseif x>xi(end-1)
L1=(x-xi(end-1))*(x-xi(end))/(xi(end-2)-xi(end-1))/(xi(end-2)-xi(end));
L2=(x-xi(end-2))*(x-xi(end))/(xi(end-1)-xi(end-2))/(xi(end-1)-xi(end));
L3=(x-xi(end-2))*(x-xi(end-1))/(xi(end)-xi(end-2))/(xi(end)-xi(end-1));
y=L1*yi(1)+L2*yi(2)+L3*yi(3);
else
for k=2:n-1
if xi(k+1)>x
if abs(x-xi(k))<abs(x-xi(k+1))
i=k-1;
else
i=k;
end
L1=(x-xi(i+1))*(x-xi(i+2))/(xi(i)-xi(i+1))/(xi(i)-xi(i+2));
L2=(x-xi(i))*(x-xi(i+2))/(xi(i+1)-xi(i))/(xi(i+1)-xi(i+2));
L3=(x-xi(i))*(x-xi(i+1))/(xi(i+2)-xi(i))/(xi(i+2)-xi(i+1));
y=L1*yi(i)+L2*yi(i+1)+L3*yi(i+2);
end
end
end
end