三次样条插值的Matlab

function [m,H,lambda,mu,D,A,dY,s]=ClampedSp_1 (X,Y,dyo,dyn)
m=length(X);A=zeros(m,m); n=m-1;H=zeros(1,n);lambda=zeros(1,n);
mu=zeros(1,n);lambda(1)=1;A(1,1)=2;A(1,2)=lambda(1);
D=zeros(1,n);mu(n)=1;
for k=1:n
    hk=X(k+1)-X(k);H(k)=hk;
end
for k=1:n-1
    lambdak=H(k+1)/(H(k)+H(k+1)); lambda(k+1)=lambdak;
    muk=1-lambda(k+1);mu(k)=muk;
    dk=6*((Y(k+2)-Y(k+1))/H(k+1)-(Y(k+1)-Y(k))/H(k))/(H(k+1)+H(k));D(k+1)=dk;
end
D(1)=6*((Y(2)-Y(1))/H(1)-dyo)/H(1);
D(m)=6*(dyn-(Y(5)-Y(4))/H(4))/H(4);
for i=1:n
    A(i,i)=2;A(m,m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i);
end
dY=A\D';
syms x 
syms S1 S2 S3 S4 
S=[S1 S2 S3 S4];
for k=1:n
    
    S(k)=dY(k)*((X(k+1)-x)^3)/6/H(k)+dY(k+1)*((x-X(k))^3)/6/H(k)+(Y(k)-dY(k)/6*H(k)^2)*(X(k+1)-x)/H(k)+(Y(k+1)-dY(k+1)/6*H(k)^2)*(x-X(k))/H(k);
    s=S';
end
end

运行结果：

>> X=[0.25 0.30 0.39 0.45 0.53];

>> Y=[0.5 0.5477 0.6245 0.6708 0.7280];
>> dyo=1.0000;

>> dyn=0.6868;

>>[m,H,dda,mu,D,A,dY,s]=JQW(X,Y,dyo,dyn)

n =

     5


H =

    0.0500    0.0900    0.0600    0.0800


dda =

    1.0000    0.6429    0.4000    0.5714


mu =

    0.3571    0.6000    0.4286    1.0000


D =

   -5.5200   -4.3143   -3.2667   -2.4286   -2.1150


A =

    2.0000    1.0000         0         0         0
    0.3571    2.0000    0.6429         0         0
         0    0.6000    2.0000    0.4000         0
         0         0    0.4286    2.0000    0.5714
         0         0         0    1.0000    2.0000


dY =

   -2.0286
   -1.4627
   -1.0333
   -0.8058
   -0.6546

 
s =
 
           (2137598125399145*conj(x))/2251799813685248 - (32937999086543005*(conj(x) - 1/4)^3)/6755399441055744 + (7613445885740335*(conj(x) - 3/10)^3)/1125899906842624 + 2373614947857091/9007199254740992
  (4290822296535025*conj(x))/5066549580791808 - (116344298039144075*(conj(x) - 3/10)^3)/60798594969501696 + (164689995432715025*(conj(x) - 39/100)^3)/60798594969501696 + 3993886489764541/13510798882111488
    (1732515684734625*conj(x))/2251799813685248 + (116344298039144075*(conj(x) - 9/20)^3)/40532396646334464 - (30242813525570275*(conj(x) - 39/100)^3)/13510798882111488 + 1463927996352807/4503599627370496
 (6421983471708925*conj(x))/9007199254740992 - (147399389720506025*(conj(x) - 9/20)^3)/108086391056891904 + (30242813525570275*(conj(x) - 53/100)^3)/18014398509481984 + 12639515432295147/36028797018963968