Matlab插值算法程序集

Matlab插值算法程序集一、插值算法1.Atkenfunctionf=Atken(x,y,x0)symst;if(length(x)==length(y))n=length(x);elsedisp('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));endy=y1;simplify(y1);endif(nargin==3)f=subs(y1(n),'t',x0);%計算插值點的函數值elsesimplify(y1(n));%化簡f=collect(y1(n));%將插值多項式展開f=vpa(f,6);%將插值多項式的系數化成6位精度的小數end2.BSamplefunctionf0=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-1if(a+i*h<=x0)&&(a+i*h+h>=x0)index=i;break;endend%找到x0所在區(qū)間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:nb(i,1)=6*y(i);endb(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.DCSfunctionf=DCS(x,y,x0)symst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數不相等！');return;endc(1)=y(1);for(i=1:n-1)for(j=i+1:n)y1(j)=(x(j)-x(i))/(y(j)-y(i));endc(i+1)=y1(i+1);y=y1;endf=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);elsef=vpa(f,6);endendend;4.DHfunctionfz=DH(x,y,x0,y0,zx,zy,zxy)n=length(x);m=length(y);fori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index_x=i;break;endend%找到x0所在區(qū)間fori=1:mif(y(i)<=y0)&&(y(i+1)>=y0)index_y=i;break;endend%找到y(tǒng)0所在區(qū)間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.DLfunctionfz=DL(x,y,Z,x0,y0,eps)n=length(x);m=length(y);fori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index_x=i;break;endend%找到x0所在區(qū)間fori=1:mif(y(i)<=y0)&&(y(i+1)>=y0)index_y=i;break;endend%找到y(tǒng)0所在區(qū)間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.DTLfunctionfz=DTL(x,y,Z,x0,y0)symsst;f=0.0;n=length(x);m=length(y);fori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index_x=i;break;endend%找到x0所在區(qū)間fori=1:mif(y(i)<=y0)&&(y(i+1)>=y0)index_y=i;break;endend%找到y(tǒng)0所在區(qū)間ifindex_x==1cx(1:3)=index_x:(index_x+2);elseifindex_x==n-1cx(1:3)=(index_x-1):(index_x+1);elseifabs(x(index_x-1)-x0)>abs(x(index_x+2)-x0)cx(1:3)=(index_x):(index_x+2);elsecx(1:3)=(index_x-1):(index_x+1);endendend%找到離x0最近的三個x坐標ifindex_y==1cy(1:3)=index_y:(index_y+2);elseifindex_y==m-1cy(1:3)=(index_y-1):(index_y+1);elseifabs(y(index_y-1)-y0)>=abs(y(index_y+2)-y0)cy(1:3)=(index_y):(index_y+2);elsecy(1:3)=(index_y-1):(index_y+1);endendend%找到離y0最近的三個y坐標fori=1:3i1=mod(i+1,3);if(i1==0)i1=3;endi2=mod(i+2,3);if(i2==0)i2=3;endforj=1:3j1=mod(j+1,3);if(j1==0)j1=3;endj2=mod(j+2,3);if(j2==0)j2=3;endf=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)));%插值多項式endendfz=subs(f,'[ts]',[x0y0]);FCZfunctionfz=DTL(x,y,Z,x0,y0)symsst;f=0.0;n=length(x);m=length(y);fori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index_x=i;break;endend%找到x0所在區(qū)間fori=1:mif(y(i)<=y0)&&(y(i+1)>=y0)index_y=i;break;endend%找到y(tǒng)0所在區(qū)間ifindex_x==1cx(1:3)=index_x:(index_x+2);elseifindex_x==n-1cx(1:3)=(index_x-1):(index_x+1);elseifabs(x(index_x-1)-x0)>abs(x(index_x+2)-x0)cx(1:3)=(index_x):(index_x+2);elsecx(1:3)=(index_x-1):(index_x+1);endendend%找到離x0最近的三個x坐標ifindex_y==1cy(1:3)=index_y:(index_y+2);elseifindex_y==m-1cy(1:3)=(index_y-1):(index_y+1);elseifabs(y(index_y-1)-y0)>=abs(y(index_y+2)-y0)cy(1:3)=(index_y):(index_y+2);elsecy(1:3)=(index_y-1):(index_y+1);endendend%找到離y0最近的三個y坐標fori=1:3i1=mod(i+1,3);if(i1==0)i1=3;endi2=mod(i+2,3);if(i2==0)i2=3;endforj=1:3j1=mod(j+1,3);if(j1==0)j1=3;endj2=mod(j+2,3);if(j2==0)j2=3;endf=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)));%插值多項式endendfz=subs(f,'[ts]',[x0y0]);Gaussfunctionf=Gauss(x,y,x0)if(length(x)==length(y))n=length(x);elsedisp('x和y的維數不相等！');return;endxx=linspace(x(1),x(n),(x(2)-x(1)));if(xx~=x)disp('節(jié)點之間不是等距的！');return;endif(mod(n,2)==1)if(nargin==2)f=GStirling(x,y,n);elseif(nargin==3)f=GStirling(x,y,n,x0);endendelseif(nargin==2)f=GBessel(x,y,n);elseif(nargin==3)f=GBessel(x,y,n,x0);endendendfunctionf=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);endif(mod(i,2)==1)c(i)=(y1((i+n)/2)+y1((i+n+2)/2))/2;elsec(i)=y1((i+n+1)/2)/2;endif(mod(i,2)==1)l=t+(i-1)/2;for(k=1:i-1)l=l*(t+(i-1)/2-k);endelsel_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);endl=l_1+l_2;endl=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)));endendendfunctionf=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);endif(mod(i,2)==1)c(i)=y1((i+n+1)/2)/2;elsec(i)=(y1((i+n)/2)+y1((i+n+2)/2))/2;endif(mod(i,2)==0)l=t+i/2-1;for(k=1:i-1)l=l*(t+i/2-1-k);endelsel_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);endl=l_1+l_2;endl=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)));endendEnd9.Hermitefunctionf=Hermite(x,y,y_1,x0)symst;f=0.0;if(length(x)==length(y))if(length(y)==length(y_1))n=length(x);elsedisp('y和y的導數的維數不相等！');return;endelsedisp('x和y的維數不相等！');return;endfori=1:nh=1.0;a=0.0;forj=1:nif(j~=i)h=h*(t-x(j))^2/((x(i)-x(j))^2);a=a+1/(x(i)-x(j));endendf=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);elsef=vpa(f,6);endendEnd10.Languagefunctionf=Language(x,y,x0)symst;if(length(x)==length(y))n=length(x);elsedisp('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);%計算插值點的函數值elsef=collect(f);%將插值多項式展開f=vpa(f,6);%將插值多項式的系數化成6位精度的小數endendend11.Nevillefunctionf=Neville(x,y,x0)symst;if(length(x)==length(y))n=length(x);elsedisp('x和y的維數不相等！');return;endy1(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));elsey1(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)));endendy=y1;if(i==n-1)if(nargin==3)f=subs(y(n-1),'t',x0);elsef=vpa(y(n-1),6);endendend12.Newtonfunctionf=Newton(x,y,x0)symst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數不相等！');return;endf=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));endc(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);elsef=collect(f);%將插值多項式展開f=vpa(f,6);endendend13.Newtonbackfunctionf=Newtonback(x,y,x0)symst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數不相等！');return;endf=y(n);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1)));if(xx~=x)disp('節(jié)點之間不是等距的！');return;endfor(i=1:n-1)for(j=i+1:n)y1(j)=y(j)-y(j-1);endc(i)=y1(n);l=t;for(k=1:i-1)l=l*(t+k);end;f=f+c(i)*l/factorial(i);simplify(f);y=y1;if(i==n-1)if(nargin==3)f=subs(f,'t',(x0-x(n))/(x(2)-x(1)));elsef=collect(f);f=vpa(f,6);endendend14.Newtonforwardfunctionf=Newtonforward(x,y,x0)symst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數不相等！');return;endf=y(1);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1)));if(xx~=x)disp('節(jié)點之間不是等距的！');return;endfor(i=1:n-1)for(j=1:n-i)y1(j)=y(j+1)-y(j);endc(i)=y1(1);l=t;for(k=1:i-1)l=l*(t-k);end;f=f+c(i)*l/factorial(i);simplify(f);y=y1;if(i==n-1)if(nargin==3)f=subs(f,'t',(x0-x(1))/(x(2)-x(1)));elsef=collect(f);f=vpa(f,6);endendend15.SecSamplefunction[f,f0,fd0]=SecSample(x,y,y_1,x0)symst;f=0.0;f0=0.0;if(length(x)==length(y))if(length(y)==length(y_1))n=length(x);elsedisp('y和y的導數的維數不相等！');return;endelsedisp('x和y的維數不相等！');return;end%維數檢查fori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endend%找到x0所在區(qū)間d=y_1(1)*(x(2)-x(1))/2+y(1);fori=2:n-1d=d+y_1(i)*(x(i+1)-x(i-1))/2;endh=x(index+1)-x(index);%x0所在區(qū)間長度f=y_1(index+1)*(t-x(index))^2/2/h+...y_1(index)*(t-x(index+1))^2/2/h+d;%x0所在區(qū)間的插值函數fd=(t-x(index))*y_1(index+1)/h+y_1(index)*(x(index+1)-t)/h;%x0所在區(qū)間的插值函數的導數f0=subs(f,'t',x0);%x0處的插值fd0=subs(fd,'t',x0);%x0處的導數插值16.SubHermitefunction[f,f0]=SubHermite(x,y,y_1,x0)symst;f=0.0;f0=0.0;if(length(x)==length(y))if(length(y)==length(y_1))n=length(x);elsedisp('y和y的導數的維數不相等！');return;endelsedisp('x和y的維數不相等！');return;end%維數檢查fori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endend%找到x0所在區(qū)間h=x(index+1)-x(index);%x0所在區(qū)間長度fl=y(index)*(1+2*(t-x(index))/h)*(t-x(index+1))^2/h/h+...y(index+1)*(1-2*(t-x(index+1))/h)*(t-x(index))^2/h/h;fr=y_1(index)*(t-x(index))*(t-x(index+1))^2/h/h+...y_1(index+1)*(t-x(index+1))*(t-x(index))^2/h/h;f=fl+fr;%x0所在區(qū)間的插值函數f0=subs(f,'t',x0);%x0處的插值17.ThrSample1function[f,f0]=ThrSample1(x,y,y_1,y_N,x0)symst;f=0.0;f0=0.0;if(length(x)==length(y))n=length(x);elsedisp('x和y的維數不相等！');return;end%維數檢查fori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endend%找到x0所在區(qū)間fori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endend%找到x0所在區(qū)間A=diag(2*ones(1,n));%求解m的系數矩陣u=zeros(n-2,1);lamda=zeros(n-1,1);c=zeros(n,1);fori=2:n-1u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+...3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));A(i,i+1)=u(i-1);A(i,i-1)=lamda(i);%形成系數矩陣及向量cendc(1)=2*y_1;c(n)=2*y_N;m=followup(A,c);%用追趕法求解方程組h=x(index+1)-x(index);%x0所在區(qū)間長度f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+...y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+...m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-...m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h;%x0所在區(qū)間的插值函數f0=subs(f,'t',x0);%x0處的插值18.ThrSample2function[f,f0]=ThrSample2(x,y,y2_1,y2_N,x0)symst;f=0.0;f0=0.0;if(length(x)==length(y))n=length(x);elsedisp('x和y的維數不相等！');return;end%維數檢查fori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endend%找到x0所在區(qū)間fori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endend%找到x0所在區(qū)間A=diag(2*ones(1,n));%求解m的系數矩陣A(1,2)=1;A(n,n-1)=1;u=zeros(n-2,1);lamda=zeros(n-1,1);c=zeros(n,1);fori=2:n-1u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+...3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));A(i,i+1)=u(i-1);A(i,i-1)=lamda(i);%形成系數矩陣及向量cendc(1)=3*(y(2)-y(1))/(x(2)-x(1))-(x(2)-x(1))*y2_1/2;c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))-(x(n)-x(n-1))*y2_N/2;m=followup(A,c