首先:
几种常用的插值方法主要有:Lagrange插值、Newton插值、Hermite插值;分段插值方法主要有:分段线性插值、分段三次Hermite插值、三次样条插值。
接下来:
已知的插值公式:
已知的分段插值公式:
将以上插值公式使用Matlab算法实现:
X为x的行向量,Y为y的行向量,Z为y的一阶导数向量或者边界一阶导数行向量。不同的方法用到的部分代码可能相同。
代码:
Lagrange插值公式求解及作图代码:
function [L]=Lagrange(X,Y)
[C,L,L1,l]=Lagran(X,Y);
X1=1:0.01:5;
Y1=polyval(C,X1);
plot(X,Y,'*',X1,Y1,'r');
grid on;
xlabel('x');
ylabel('y');
function [C,L,L1,l]=Lagran(X,Y)
m=length(X);
L=ones(m,m);
for k=1:m
v=1;
for i=1:m
if k~=i
v=conv(v,poly(X(i)))/(X(k)-X(i));
end
end
L1(k,:)=v;
l(k,:)=poly2sym(v);
end
C=Y*L1;
L=Y*l;
Newton插值公式求解及作图代码:
function [L]=Newton(X,Y)
[A,C,L]=Newploy(X,Y);
X1=1:0.01:5;
Y1=polyval(C,X1);
plot(X,Y,'*',X1,Y1,'r');
grid on
xlabel('x');
ylabel('y');
function [A,C,L]=Newploy(X,Y)
n=length(X);
A=zeros(n,n);
A(:,1)=Y';
s=0.0;
p=1.0;
q=1.0;
c1=1.0;
for j=2:n
for i=j:n
A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));
end
b=poly(X(j-1));
q1=conv(q,b);
c1=c1*j;
q=q1;
end
C=A(n,n);
b=poly(X(n));
q1=conv(q1,b);
for k=(n-1):-1:1
C=conv(C,poly(X(k)));
d=length(C);
C(d)=C(d)+A(k,k);
end
L(k,:)=poly2sym(C);
Hermite插值公式求解及作图代码:
function [L]=Hermite(X,Y,Z)
[A,C,L]=Herm1(X,Y,Z);
X1=1:0.01:5;
Y1=polyval(C,X1);
plot(X,Y,'*',X1,Y1,'r');
grid on
xlabel('x');
ylabel('y');
function [A,C,L]=Herm1(X,Y,Z)
s=length(X);
n=2*s;
X2=zeros(1,n);
Y2=zeros(1,n);
Z2=zeros(1,n);
A=zeros(n,n);
for m=1:s
X2(2*m-1)=X(m);
X2(2*m)=X(m);
Y2(2*m-1)=Y(m);
Y2(2*m)=Y(m);
Z2(2*m)=Z(m);
end
A(:,1)=Y2';
s=0.0;
p=1.0;
q=1.0;
c1=1.0;
for j=2:n
for i=j:n
if X2(i)==X2(i-j+1)
A(i,j)=Z2(i);
else
A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X2(i)-X2(i-j+1));
end
end
b=poly(X2(j-1));
q1=conv(q,b);
c1=c1*j;
q=q1;
end
C=A(n,n);
b=poly(X2(n));
q1=conv(q1,b);
for k=(n-1):-1:1
C=conv(C,poly(X2(k)));
d=length(C);
C(d)=C(d)+A(k,k);
end
L(k,:)=poly2sym(C);
分段线性插值公式求解及作图代码:
function FD(X,Y)
n=length(X);
X1=X;
Y1=Y;
for i=1:n-1
x=X1(i):0.01:X1(i+1);
y=Y1(i)*(x-X1(i+1))/(X1(i)-X1(i+1))+Y1(i+1)*(x-X1(i))/(X1(i+1)-X1(i));
plot(X,Y,'*',X1,Y1,'r');
hold on
end
xlabel('x');
ylabel('y');
syms x;
for i=1:n-1
y=Y1(i)*(x-X1(i+1))/(X1(i)-X1(i+1))+Y1(i+1)*(x-X1(i))/(X1(i+1)-X1(i))
end
分段三次Hermite插值公式求解及作图代码:
function TH(X,Y,Z)
n=length(X);
X3=zeros(2);
Y3=zeros(2);
Z3=zeros(2);
for i=1:n-1
X3(1)=X(i);
X3(2)=X(i+1);
Y3(1)=Y(i);
Y3(2)=Y(i+1);
Z3(1)=Z(i);
Z3(2)=Z(i+1);
[A,C,L]=Herm1(X3,Y3,Z3)
X1=X3(1):0.01:X3(2);
Y1=polyval(C,X1);
plot(X,Y,'*',X1,Y1,'r');
grid on
hold on
end
xlabel('x');
ylabel('y');
三次样条插值公式求解及作图代码:
function SY(X,Y,Z)
n=length(X);
Y1=zeros(n,1);
A=2*eye(n);
for i=2:n-1
A(i,i-1)=(X(i)-X(i-1))/(X(i+1)-X(i-1));
A(i,i+1)=1-(X(i)-X(i-1))/(X(i+1)-X(i-1));
end
A(1,2)=1;
A(n,n-1)=1;
Y1(1)=6/(X(2)-X(1))*((Y(2)-Y(1))/(X(2)-X(1))-Z(1));
Y1(n)=6/(X(n)-X(n-1))*(Z(2)-(Y(n)-Y(n-1))/(X(n)-X(n-1)));
for i=2:n-1
Y1(i)=6/(X(i+1)-X(i-1))*((Y(i+1)-Y(i))/(X(i+1)-X(i))-(Y(i)-Y(i-1))/(X(i)-X(i-1)));
end
U=A;
L=eye(n);
for k=1:n-1
L(k+1:n,k)=U(k+1:n,k)/U(k,k);
U(k+1:n,k+1:n)=U(k+1:n,k+1:n)-L(k+1:n,k)*U(k,k+1:n);
U(k+1:n,k)=zeros(n-k,1);
end
for j=1:n-1
Y1(j)=Y1(j)/L(j,j);
Y1(j+1:n)=Y1(j+1:n)-Y1(j)*L(j+1:n,j);
end
Y1(n)=Y1(n)/L(n,n);
for j=n:-1:2
Y1(j)=Y1(j)/U(j,j);
Y1(1:j-1)=Y1(1:j-1)-Y1(j)*U(1:j-1,j);
end
Y1(1)=Y1(1)/U(1,1);
for i=1:n-1
x=X(i):0.01:X(i+1);
y=Y1(i)/6*(X(i+1)-x).^3/(X(i+1)-X(i))+Y1(i+1)/6*(x-X(i)).^3/(X(i+1)-X(i))+(Y(i)-Y1(i)/6*(X(i+1)-X(i)).^2)*(X(i+1)-x)/(X(i+1)-X(i))+(Y(i+1)-Y1(i+1)/6*(X(i+1)-X(i)).^2)*(x-X(i))/(X(i+1)-X(i));
plot(X,Y,'*',x,y,'r');
hold on
end
xlabel('x');
ylabel('y');
syms x;
for i=1:n-1
y=Y1(i)/6*(X(i+1)-x).^3/(X(i+1)-X(i))+Y1(i+1)/6*(x-X(i)).^3/(X(i+1)-X(i))+(Y(i)-Y1(i)/6*(X(i+1)-X(i)).^2)*(X(i+1)-x)/(X(i+1)-X(i))+(Y(i+1)-Y1(i+1)/6*(X(i+1)-X(i)).^2)*(x-X(i))/(X(i+1)-X(i))
end
对同一例题采用以上方法分别求解:
Lagrange插值:
Newton插值:
Hermite插值:
分段线性插值:
分段三次Hermite插值:
三次样条插值:
解释:
Lagrange插值和Hermite插值所得结果相同的原因是:两者均采用n次多项式插值,同属于代数插值的范畴。