Matlab插值算法程序集.docx

1、Matlab插值算法程序集 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);endif(nargin = 3)

2、 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) index = i; break; endend %找到x0所在区间A = diag(4*

3、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);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-

4、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

5、 = DCS(x,y,x0)syms t;if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp(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); end c(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 =

6、subs(f,t,x0); else f = vpa(f,6); end endend;4.DHfunction fz = DH(x,y,x0,y0,zx,zy,zxy)n = length(x);m = length(y);for i=1:n if(x(i)=x0) index_x = i; break; endend %找到x0所在区间for i=1:m if(y(i)=y0) index_y = i; break; endend %找到y0所在区间hx = x(index_x+1) - x(index_x);hy = y(index_y+1) - y(index_y);tx = (x0

7、- 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, in

8、dex_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.DLfunction

9、 fz = DL(x,y,Z,x0,y0,eps)n = length(x);m = length(y);for i=1:n if(x(i)=x0) index_x = i; break; endend %找到x0所在区间for i=1:m if(y(i)=y0) index_y = i; break; endend %找到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+

10、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.DTLfunction fz = DTL(x,y,Z,x0,y0)syms s t;f = 0.0;n = length(x);m = length(y);for i=1:n if(

11、x(i)=x0) index_x = i; break; endend %找到x0所在区间for i=1:m if(y(i)=y0) index_y = i; break; endend %找到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)

12、= (index_x-1):(index_x+1); end endend %找到离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 endend %找到离y0最近的三个y

13、坐标for i=1:3 i1 = mod(i+1,3); if(i1 = 0) i1 = 3; end i2 = mod(i+2,3); if(i2 = 0) i2 = 3; end for j=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(

14、cy(j1)/(y(cy(j)-y(cy(j2); %插值多项式 endendfz = subs(f,t s,x0 y0);7.FCZfunction 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) index_x = i; break; endend %找到x0所在区间for i=1:m if(y(i)=y0) index_y = i; break; endend %找到y0所在区间if index_x = 1 cx(1:3) = index_x:(index_x+

15、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 endend %找到离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):(ind

16、ex_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 endend %找到离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) i2 = 3; end for j=1:3 j1 = mod(j+1,3); if(j1 = 0) j1 = 3; end j2

17、= 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); %插值多项式 endendfz = subs(f,t s,x0 y0);8.Gaussfunction f = Gauss(x,y,x0)if(length(x) = length(y) n = length(x);else disp(x

18、和y的维数不相等！); return;endxx =linspace(x(1),x(n),(x(2)-x(1);if(xx = x) disp(节点之间不是等距的！); return;endif( mod(n,2) =1) if(nargin = 2) f = GStirling(x,y,n); else if(nargin = 3) f = GStirling(x,y,n,x0); end endelse if(nargin = 2) f = GBessel(x,y,n); else if(nargin = 3) f = GBessel(x,y,n,x0); end endendfuncti

19、on f = GStirling(x,y,n,x0)syms t;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; 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

20、) 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 endendfunction f = GBessel(x,y,n,x0)syms t;nn = n/2;f = (y(nn)+y(nn+1)/2;for(i=1:n-1) for(j=i+1

21、: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;

22、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 endEnd9.Hermitefunction f = Hermite(x,y,y_1,x0)syms t;f = 0.0;if(length(x) = length(y) if(length(y) = length(y_1) n = length(x); else disp(y和y的导数的维数不相等！); retur

23、n; endelse disp(x和y的维数不相等！); return;endfor i=1:n h = 1.0; a = 0.0; for j=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 endEnd10.Languagefunction f = Language(x,y,x0)s

24、yms t;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); %计

25、算插值点的函数值 else f = collect(f); %将插值多项式展开 f = vpa(f,6); %将插值多项式的系数化成6位精度的小数 end endend11.Nevillefunction f = Neville(x,y,x0)syms t;if(length(x) = length(y) n = length(x);else disp(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)-

26、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 end12.Newtonfunction f = Newton(x,y,x0)syms t;if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp(x和y

27、的维数不相等！); 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); 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 endend13.Newtonbackfunction f = New

28、tonback(x,y,x0)syms t;if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp(x和y的维数不相等！); return;endf = y(n);y1 = 0;xx =linspace(x(1),x(n),(x(2)-x(1);if(xx = x) disp(节点之间不是等距的！); return;endfor(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