三次样条插值.doc

《三次样条插值.doc》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、三次样条插值2010-05-2810:14functionvarargout=saxplaxliu(varargin)clear,clcx=[0.250.30.390.450.53];y=[0.50.54770.62540.67080.728];fd=[10.6868];x0=0.4;yx0=splineL(x,y,fd,x0)functionyx0=splineL(x,y,fd,x0)%x=x0x1xn%y=y0y1yn%fd=[df0dfn];%n=length(x)-1;f=x;h=x;m=x;fori=1:nh(i)=

2、x(i+1)-x(i);f(i)=(y(i+1)-y(i))/h(i);endB=2*ones(1,n-1);m(1)=fd(1);m(end)=fd(2);A=B;C=B;g=A;fori=1:n-1C(i)=h(i)/(h(i+1)+h(i));A(i)=h(i+1)/(h(i+1)+h(i));g(i)=3*(A(i)*f(i)+C(i)*f(i+1));endfL(1)=g(1)-A(1)*fd(1);fL(n-1)=g(n-1)-C(n-1)*fd(2);fL(2:n-2)=g(2:n-2);m(2:n)=zgf(A

3、,B,C,fL);yx0=fliu(x,y,m,x0);functiony=fliu(xj,yj,mj,x)n=length(xj);y=0;fori=1:ns=1;forj=1:nifi~=js=s*(x-xj(j))/(xj(i)-xj(j));endendbi=(x-xj(i))*s^2;s1=0;forj=1:nifi~=js1=s1+1/(xj(i)-xj(j));endendy=y+yj(i)*(1-2*(x-xj(i))*s1)*s^2+mj(i)*bi;endfunctionx=zgf(A,B,C,f)%解[b

4、1c1%   a2b2c2%      ..  .       *x=f%                 %              anbn]n=length(B);B1=zeros(1,n-1);Y=zeros(1,n);x1=zeros(1,n);B1(1)=C(1)/B(1);fori=2:n-1B1(i)=C(i)/(B(i)-A(i)*B1(i-1));endY(1)=f(1)/B(1);fori=2:nY(i)=(f(i)-A(i)*Y(i-1))/(B(i)-A(i)*B1(i-1));endx1(n)=Y

5、(n);fori=n-1:-1:1x1(i)=Y(i)-B1(i)*x1(i+1);   endx=x1;