Matlab插值算法程序集.docx

Matlab插值算法程序集.docx

《Matlab插值算法程序集.docx》由会员分享，可在线阅读，更多相关《Matlab插值算法程序集.docx（27页珍藏版）》请在冰豆网上搜索。

Matlab插值算法程序集

Matlab插值算法程序集

一、插值算法

1.Atken

functionf=Atken（x,y,x0）

symst;

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

functionf0=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）;

fori=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;%形成系数矩阵

fori=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

functionf=DCS（x,y,x0）

symst;

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

functionfz=DH（x,y,x0,y0,zx,zy,zxy）

n=length（x）;

m=length（y）;

fori=1:

n

if（x（i）<=x0）&&（x（i+1）>=x0）

index_x=i;

break;

end

end%找到x0所在区间

fori=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）^2tx*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

functionfz=DL（x,y,Z,x0,y0,eps）

n=length（x）;

m=length（y）;

fori=1:

n

if（x（i）<=x0）&&（x（i+1）>=x0）

index_x=i;

break;

end

end%找到x0所在区间

fori=1:

m

if（y（i）<=y0）&&（y（i+1）>=y0）

index_y=i;

break;

end

end%找到y0所在区间

A=[1x（index_x）y（index_y）x（index_x）*y（index_y）;

1x（index_x+1）y（index_y+1）x（index_x+1）*y（index_y+1）;

1x（index_x）y（index_y+1）x（index_x）*y（index_y+1）;

1x（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=[1x0y0x0*y0]*B;

6.DTL

functionfz=DTL（x,y,Z,x0,y0）

symsst;

f=0.0;

n=length（x）;

m=length（y）;

fori=1:

n

if（x（i）<=x0）&&（x（i+1）>=x0）

index_x=i;

break;

end

end%找到x0所在区间

fori=1:

m

if（y（i）<=y0）&&（y（i+1）>=y0）

index_y=i;

break;

end

end%找到y0所在区间

ifindex_x==1

cx（1:

3）=index_x:

（index_x+2）;

else

ifindex_x==n-1

cx（1:

3）=（index_x-1）:

（index_x+1）;

else

ifabs（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坐标

ifindex_y==1

cy（1:

3）=index_y:

（index_y+2）;

else

ifindex_y==m-1

cy（1:

3）=（index_y-1）:

（index_y+1）;

else

ifabs（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坐标

fori=1:

3

i1=mod（i+1,3）;

if（i1==0）

i1=3;

end

i2=mod（i+2,3）;

if（i2==0）

i2=3;

end

forj=1:

3

j1=mod（j+1,3）;

if（j1==0）

j1=3;

end

j2=mod（j+2,3）;

if（j2==0）

j2=3;

end

f=f+Z（cx（i）,cy（j））*（（t-x（cx（i1）））*（t-x（cx（i2）））/（x（cx（i））-x（cx（i1）））/（x（cx（i））-x（cx（i2））））*...

（s-y（cy（j1）））*（s-y（cy（j2）））/（y（cy（j））-y（cy（j1）））/（y（cy（j））-y（cy（j2）））;

%插值多项式

end

end

fz=subs（f,'[ts]',[x0y0]）;

7.FCZ

functionfz=DTL（x,y,Z,x0,y0）

symsst;

f=0.0;

n=length（x）;

m=length（y）;

fori=1:

n

if（x（i）<=x0）&&（x（i+1）>=x0）

index_x=i;

break;

end

end%找到x0所在区间

fori=1:

m

if（y（i）<=y0）&&（y（i+1）>=y0）

index_y=i;

break;

end

end%找到y0所在区间

ifindex_x==1

cx（1:

3）=index_x:

（index_x+2）;

else

ifindex_x==n-1

cx（1:

3）=（index_x-1）:

（index_x+1）;

else

ifabs（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坐标

ifindex_y==1

cy（1:

3）=index_y:

（index_y+2）;

else

ifindex_y==m-1

cy（1:

3）=（index_y-1）:

（index_y+1）;

else

ifabs（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坐标

fori=1:

3

i1=mod（i+1,3）;

if（i1==0）

i1=3;

end

i2=mod（i+2,3）;

if（i2==0）

i2=3;

end

forj=1:

3

j1=mod（j+1,3）;

if（j1==0）

j1=3;

end

j2=mod（j+2,3）;

if（j2==0）

j2=3;

end

f=f+Z（cx（i）,cy（j））*（（t-x（cx（i1）））*（t-x（cx（i2）））/（x（cx（i））-x（cx（i1）））/（x（cx（i））-x（cx（i2））））*...

（s-y（cy（j1）））*（s-y（cy（j2）））/（y（cy（j））-y（cy（j1）））/（y（cy（j））-y（cy（j2）））;

%插值多项式

end

end

fz=subs（f,'[ts]',[x0y0]）;

8.Gauss

functionf=Gauss（x,y,x0）

if（length（x）==length（y））

n=length（x）;

else

disp（'x和y的维数不相等！

'）;

return;

end

xx=linspace（x

（1）,x（n）,（x

（2）-x

（1）））;

if（xx~=x）

disp（'节点之间不是等距的！

'）;

return;

end

if（mod（n,2）==1）

if（nargin==2）

f=GStirling（x,y,n）;

elseif（nargin==3）

f=GStirling（x,y,n,x0）;

end

end

else

if（nargin==2）

f=GBessel（x,y,n）;

elseif（nargin==3）

f=GBessel（x,y,n,x0）;

end

end

end

functionf=GStirling（x,y,n,x0）

symst;

nn=（n+1）/2;

f=y（nn）;

for（i=1:

n-1）

for（j=i+1:

n）

y1（j）=y（j）-y（j-1）;

end

if（mod（i,2）==1）

c（i）=（y1（（i+n）/2）+y1（（i+n+2）/2））/2;

else

c（i）=y1（（i+n+1）/2）/2;

end

if（mod（i,2）==1）

l=t+（i-1）/2;

for（k=1:

i-1）

l=l*（t+（i-1）/2-k）;

end

else

l_1=t+i/2-1;

l_2=t+i/2;

for（k=1:

i-1）

l_1=l_1*（t+i/2-1-k）;

l_2=l_2*（t+i/2-k）;

end

l=l_1+l_2;

end

l=l/factorial（i）;

f=f+c（i）*l;

simplify（f）;

f=vpa（f,6）;

y=y1;

if（i==n-1）

if（nargin==4）

f=subs（f,'t',（x0-x（nn））/（x

（2）-x

（1）））;

end

end

end

functionf=GBessel（x,y,n,x0）

symst;

nn=n/2;

f=（y（nn）+y（nn+1））/2;

for（i=1:

n-1）

for（j=i+1:

n）

y1（j）=y（j）-y（j-1）;

end

if（mod（i,2）==1）

c（i）=y1（（i+n+1）/2）/2;

else

c（i）=（y1（（i+n）/2）+y1（（i+n+2）/2））/2;

end

if（mod（i,2）==0）

l=t+i/2-1;

for（k=1:

i-1）

l=l*（t+i/2-1-k）;

end

else

l_1=t+（i-1）/2;

l_2=t+（i-1）/2-1;

for（k=1:

i-1）

l_1=l_1*（t+（i-1）/2-k）;

l_2=l_2*（t+（i-1）/2-1-k）;

end

l=l_1+l_2;

end

l=l/factorial（i）;

f=f+c（i）*l;

simplify（f）;

f=vpa（f,6）;

y=y1;

if（i==n-1）

if（nargin==4）

f=subs（f,'t',（x0-x（nn））/（x

（2）-x

（1）））;

end

end

End

9.Hermite

functionf=Hermite（x,y,y_1,x0）

symst;

f=0.0;

if（length（x）==length（y））

if（length（y）==length（y_1））

n=length（x）;

else

disp（'y和y的导数的维数不相等！

'）;

return;

end

else

disp（'x和y的维数不相等！

'）;

return;

end

fori=1:

n

h=1.0;

a=0.0;

forj=1:

n

if（j~=i）

h=h*（t-x（j））^2/（（x（i）-x（j））^2）;

a=a+1/（x（i）-x（j））;

end

end

f=f+h*（（x（i）-t）*（2*a*y（i）-y_1（i））+y（i））;

if（i==n）

if（nargin==4）

f=subs（f,'t',x0）;

else

f=vpa（f,6）;

end

end

End

10.Language

functionf=Language（x,y,x0）

symst;

if（length（x）==length（y））

n=length（x）;

else

disp（'x和y的维数不相等！

'）;

return;

end%检错

f=0.0;

for（i=1:

n）

l=y（i）;

for（j=1:

i-1）

l=l*（t-x（j））/（x（i）-x（j））;

end;

for（j=i+1:

n）

l=l*（t-x（j））/（x（i）-x（j））;%计算拉格朗日基函数

end;

f=f+l;%计算拉格朗日插值函数

simplify（f）;%化简

if（i==n）

if（nargin==3）

f=subs（f,'t',x0）;%计算插值点的函数值

else

f=collect（f）;%将插值多项式展开

f=vpa（f,6）;%将插值多项式的系数化成6位精度的小数

end

end

end

11.Neville

functionf=Neville（x,y,x0）

symst;

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）

if（j==2）

y1（j）=y（j）+（y（j）-y（j-1））/（（t-x（j-i））/（t-x（j）））*（1-（y（j）-y（j-1））/y（j））;

else

y1（j）=y（j）+（y（j）-y（j-1））/（（t-x（j-i））/（t-x（j）））*（1-（y（j）-y（j-1））/（y（j）-y（j-2）））;

end

end

y=y1;

if（i==n-1）

if（nargin==3）

f=subs（y（n-1）,'t',x0）;

else

f=vpa（y（n-1）,6）;

end

end

end

12.Newton

functionf=Newton（x,y,x0）

symst;

if（length（x）==length（y））

n=length（x）;

c（1:

n）=0.0;

else

disp（'x和y的维数不相等！

'）;

return;

end

f=y

（1）;

y1=0;

l=1;

for（i=1:

n-1）

for（j=i+1:

n）

y1（j）=（y（j）-y（i））/（x（j）-x（i））;

end

c（i）=y1（i+1）;

l=l*（t-x（i））;

f=f+c（i）*l;

simplify（f）;

y=y1;

if（i==n-1）

if（nargin==3）

f=subs（f,'t',x0）;

else

f=collect（f）;%将插值多项式展开

f=vpa（f,6）;

end

end

end

13.Newtonback

functionf=Newtonback（x,y,x0）

symst;

if（length（x）==length（y））

n=length（x）;

c（1:

n）=0.0;

else

disp（'x和y的维数不相等！

'）;

return;

end

f=y（n）;

y1=0;

xx=linspace（x

（1）,x（n）,（x

（2）-x

（1）））;

if（xx~=x）

disp（'节点之间不是等距的！

'）;

return;

end

for（i=1:

n-1）

for（j=i+1:

n）

y1（j）=y（j）-y（j-1）;

end

c（i）=y1（n）;

l=t;

for（k=1:

i-1）

l=l*（t+k）;

end;

f=f+c（i）*l/factorial（i）;

simplify（f