B样条基函数求导

clc;clear;close

p=3;
NodeVector=[0,0,0,0,1,2,4,7,7,7,7];
u=4.5;
i=7;
der=1;//求导次数

[B,D]= MyBaseFun(i, p , u, NodeVector);

MyDersBasisFuns( i,u,p,B,D,der,NodeVector )

 

 

function [ Ders ] = MyDersBasisFuns( i,u,p,B,D,der,U )
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
for dp=1:p+1
    Ders(1,i-p-1+dp)=B(i-p-1+dp);
end

for r=1:p+1%%这个表示求个B的个数
  a(1,1)=1;
  jp=i-p+r-1;
  for k=2:der+1%%k表示从2开始求一阶导
    a(k,1)=a(k-1,1)/(U(jp+p-k+2)-U(jp));
    for d=2:k-1%%里面求aki用的
        if (U(jp+p+d-1-k+2)-U(jp+d-1))==0
            a(k,d)=0;
        else
            a(k,d)=(a(k-1,d)-a(k-1,d-1))/(U(jp+p+d-1-k+2)-U(jp+d-1));
        end
        
    end
    
    if (U(jp+p+1)-U(jp+k-1))==0
        a(k,k)=0;
    else
        a(k,k)=(-a(k-1,k-1))/(U(jp+p+1)-U(jp+k-1));
    end
    
      p1=prod(1:p);
      p2=prod(1:p-k+1);
      p3=p1/p2;
      p4=0;
    for f=1:k
        if jp+f-1>size(D,1)
            D(jp+f-1,p-k+2)=0;
        end
       p4=p4+a(k,f)*D(jp+f-1,p-k+2) 
    end
    Ders(k,jp)=p3*p4;
  end
  
end

end