Matlab 版 (精华区)

发信人: zjliu (秋天的萝卜), 信区: Matlab
标  题: 关于B-样条的一点算法
发信站: BBS 哈工大紫丁香站 (Tue May 25 16:00:27 2004)

转于饮水思源站

二次B-样条函数：

function [c, h]= ipbspline2(t,y)
%Return the coefficients of the quadratic spline interpolating
%polynomial s(x)with knots t and values y.
%And return step size h
n=length(t); %t0, t1, ..., tn
for j=2:n
   h(j)=t(j)-t(j-1);
end ;
h(1)=h(2);
h(n+1)=h(n);
del=-1; gma=2*y(1) ;
p=del*gma; q = 2;
for j = 2:n-1
   r=h(j+1)/h(j);
   del = -r*del ;
   gma = -r*gma+(r+1)*y(j);
   p = p+gma*del;
   p = p+gma*del;
   q = q+del^2 ;
end  %end j
%
c(1)=-p/q ;
for j = 2:n
   c(j) = ((h(j-1)+h(j))*y(j-1)-h(j)*c(j-1))/h(j-1);
end %end j

                 type ipbspline2x.m

function s2x= ipbspline2x(t,c,h,x)
%return the value of the quadratic spline interpolating
%polynomial s(x)with coefficients c at x
%
n=length(t); %t0, t1, ...,tn
for i= 1:length(x)
  for j = n-2:-1:1
   if x(i)-t(j) >= 0 , break, end ;
  end  %end j
%
  j=j+1;
  d=(c(j+1)*(x(i)-t(j-1))+c(j)*(t(j)-x(i)+h(j+1)))/(h(j)+h(j+1));
  d=(c(j+1)*(x(i)-t(j-1))+c(j)*(t(j)-x(i)+h(j+1)))/(h(j)+h(j+1));
  e=(c(j)*(x(i)-t(j-1)+h(j-1))+c(j-1)*(t(j-1)-x(i)+h(j)))/(h(j-1)+h(j));
  s2x(i)=(d*(x(i)-t(j-1))+e*(t(j)-x(i)))/h(j);
end %end i


eg。
t = [0 1.7 2.3 3.0 5.2 6.4 7.6 8.0];
            y = [-0.8 0.59 0.28 1.44 -1.27 -0.79 1.0 0.0] ;
            n = 7;
            x = 4.0 ;
            [c,h] = ipbspline2(t, y);
  t1=0:0.2:1.7; t2=1.7:0.1:2.3 ; t3=2.3:0.1:3.0;
  t4=3.0:0.2:5.2 ; t5=5.2:0.2:6.4 ;
  t6=6.4:0.2:7.6 ; t7=7.6:0.1:8.0 ;
       x = [t1 t2 t3 t4 t5 t6 t7] ;
            ipbspx= ipbspline2x(t,c,h, x) ;
  plot(x, ipbspx, t, y,'哦')

   http://bbs.sjtu.edu.cn/showfile?name=1069247949plot.jpg


我已经画出了f（x）=1./1+25*(x.^2)的B-spline，还可以吧，就是不是三次阿
  t6=6.4:0.2:7.6 ; t7=7.6:0.1:8.0 ;
要不干脆就LU分解那个距阵，算出来得了，呵呵！
--
╔═══════════════════╗
║★★★★★友谊第一  比赛第二★★★★★║
╚═══════════════════╝


※ 来源:·哈工大紫丁香 http://bbs.hit.edu.cn·[FROM: 202.118.229.*]