Matlab插值算法程序集 一、插值算法
1.Atken
function f = Atken(x,y,x0) syms t;
if(length(x) == length(y)) n = length(x); else
disp('x和y的维数不相等！'); return;
end %检错
y1(1:n) = t; %符号函数数组要赋初值 for(i=1:n-1) for(j=i+1:n)
y1(j) = y(j)*(t-x(i))/(x(j)-x(i))+y(i)*(t-x(j))/(x(i)-x(j)); end y = y1;
simplify(y1); end
if(nargin == 3)
f = subs(y1(n),'t',x0); %计算插值点的函数值 else
simplify(y1(n)); %化简
f = collect(y1(n)); %将插值多项式展开
f = vpa(f,6); %将插值多项式的系数化成6位精度的小数 end
2.BSample
function f0 = BSample(a,b,n,y,y_1,y_N,x0) f0 = 0.0; h = (b-a)/n;
c = zeros(n+3,1); b = zeros(n+1,1);
for i=0:n-1
if(a+i*h<=x0) && (a+i*h+h>=x0) index = i; break; end
end %找到x0所在区间
A = diag(4*ones(n+1,1)); I = eye(n+1,n+1);
AL = [I(2:n+1,:);zeros(1,n+1)]; AU = [zeros(1,n+1);I(1:n,:)];
A = A+AL+AU; %形成系数矩阵 for i=2:n
b(i,1) = 6*y(i); end
b(1) = 6*y(1)+2*h*y_1; b(n+1) = 6*y(n+1)-2*h*y_N;
d = followup(A,b); %用追赶法求出系数 c(2:n+2) = d;
c(1) = c(2) - 2*h*y_1; %c(-1) c(n+3) = c(3)+2*h*y_N; %c(n+1)
x1 = (a+index*h-h-x0)/h;
m1 = c(index+1)*(-((abs(x1))^3)/6+(x1)^2-2*abs(x1)+4/3); x2 = (a+index*h-x0)/h;
m2 = c(index+2)*((abs(x2))^3/2-(x2)^2+2/3); x3 = (a+index*h+h-x0)/h;
m3 = c(index+3)*((abs(x3))^3/2-(x3)^2+2/3); x4 = (a+index*h+2*h-x0)/h;
m4 = c(index+4)*(-((abs(x4))^3)/6+(x4)^2-2*abs(x4)+4/3); f0 = m1+m2+m3+m4; %求出插值
3.DCS
function f = DCS(x,y,x0) syms t;
if(length(x) == length(y)) n = length(x); c(1:n) = 0.0; else
disp('x和y的维数不相等！'); return; end
c(1) = y(1); for(i=1:n-1) for(j=i+1:n)
y1(j) = (x(j)-x(i))/(y(j)-y(i)); end
c(i+1) = y1(i+1);
y = y1; end
f = c(n); for(i=1:n-1)
f = c(n-i) + (t-x(n-i))/f; f = vpa(f,6); if(i==n-1)
if(nargin == 3)
f = subs(f,'t',x0); else
f = vpa(f,6); end end end;
4.DH
function fz = DH(x,y,x0,y0,zx,zy,zxy) n = length(x); m = length(y);
for i=1:n
if(x(i)<=x0)&& (x(i+1)>=x0) index_x = i; break; end
end %找到x0所在区间 for i=1:m
if(y(i)<=y0)&& (y(i+1)>=y0) index_y = i; break; end
end %找到y0所在区间 hx = x(index_x+1) - x(index_x); hy = y(index_y+1) - y(index_y);
tx = (x0 - x(index_x))/hx; %插值坐标归一化 ty = (y0 - y(index_y))/hy; %插值坐标归一化
Hl = [(1-tx)^2*(1+2*tx) tx*tx*(3-2*tx) tx*(1-tx)^2 tx*tx*(tx-1)]; %左向量 Hr = [(1-ty)^2*(1+2*ty); ty*ty*(3-2*ty); ty*(1-ty)^2 ; ty*ty*(ty-1)]; %右向量 C = [Z(index_x, index_y) Z(index_x,index_y+1) zy(index_x, index_y) ... zy(index_x, index_y+1);
Z(index_x+1, index_y) Z(index_x+1,index_y+1) zy(index_x+1, index_y) ... zy(index_x+1, index_y+1);
zx(index_x, index_y) zy(index_x, index_y+1) zxy(index_x, index_y) ... zxy(index_x, index_y+1);
zx(index_x+1, index_y) zy(index_x+1, index_y+1) zxy(index_x+1, index_y) ... zxy(index_x+1, index_y+1)]; %C矩阵 fz = Hl*C*Hr;
5.DL
function fz = DL(x,y,Z,x0,y0,eps) n = length(x); m = length(y);
for i=1:n
if(x(i)<=x0)&& (x(i+1)>=x0) index_x = i; break; end
end %找到x0所在区间 for i=1:m
if(y(i)<=y0)&& (y(i+1)>=y0) index_y = i; break; end
end %找到y0所在区间
A = [1 x(index_x) y(index_y) x(index_x)* y(index_y);
1 x(index_x+1) y(index_y+1) x(index_x+1)* y(index_y+1); 1 x(index_x) y(index_y+1) x(index_x)* y(index_y+1); 1 x(index_x+1) y(index_y) x(index_x+1)* y(index_y)]; iA = inv(A);
B = iA*[Z(index_x,index_y); Z(index_x+1,index_y+1); Z(index_x,index_y+1); Z(index_x+1,index_y)]; fz = [1 x0 y0 x0*y0]*B;
6.DTL
function fz = DTL(x,y,Z,x0,y0) syms s t; f = 0.0;
n = length(x); m = length(y);
for i=1:n
if(x(i)<=x0)&& (x(i+1)>=x0) index_x = i; break;
end
end %找到x0所在区间 for i=1:m
if(y(i)<=y0)&& (y(i+1)>=y0) index_y = i; break; end
end %找到y0所在区间
if index_x == 1
cx(1:3) = index_x:(index_x+2); else
if index_x == n-1
cx(1:3) = (index_x-1):(index_x+1); else
if abs(x(index_x-1)-x0)>abs(x(index_x+2)-x0) cx(1:3) = (index_x):(index_x+2); else
cx(1:3) = (index_x-1):(index_x+1); end end
end %找到离x0最近的三个x坐标
if index_y == 1
cy(1:3) = index_y:(index_y+2); else
if index_y == m-1
cy(1:3) = (index_y-1):(index_y+1); else
if abs(y(index_y-1)-y0)>=abs(y(index_y+2)-y0) cy(1:3) = (index_y):(index_y+2); else
cy(1:3) = (index_y-1):(index_y+1); end end
end %找到离y0最近的三个y坐标
for i=1:3 i1 = mod(i+1,3); if(i1 == 0) i1 = 3; end
i2 = mod(i+2,3); if(i2 == 0)