%姓名：汪建強 學號：1050729035 班級：B0507291
clear
A=ones(1,21);
B=ones(1,21);
X=sym('X');
J=1;F=0;
%--------------------------------------------------------------------------
for K=2:0.1:4 %---求解擬合前函數的值
    A(J)=K;
    B(J)=log(K);
    J=J+1;
end
%------------------------------------------------------------------------
for I=1:21
    C(I)=sym('1');  N=1;%－－－定義‘X’位符號變量－－－
    for J=1:21
            if(N<22)&(N~=I)   %－－－－在N不等于I的范圍內求解出系數－－
                C(I)=C(I)*(X-A(N))/(A(I)-A(N));
                N=N+1;
            else if (N==I)&(N<21)
              N=N+1;%----當I＝N時N調向下一位
            else
                    break;
                end
            end
      end
end
%-----------------------------------------------------------------------
for I=1:21  %－－－求解擬合函數的總的表達式
   F=F+C(I)*B(I);
   if I==21
       F=F
   end
end
%---------------------------------樣條擬合
U=eye(21,21);
V=ones(21,1);
Z=ones(21,1);
for I=2:20
   U(I,I-1)=(A(I)-A(I-1))/6;
   U(I,I)=(A(I+1)-A(I-1))/3;
   U(I,I+1)=(A(I+1)-A(I))/6;
   Z(I)=(B(I+1)-B(I))/(A(I+1)-A(I))-(B(I)-B(I-1))/(A(I)-A(I-1));
end
Z(1,1)=0;Z(21,1)=0;
V=inv(U)*Z;  %將U求解其逆矩陣后與R相乘。
for I=2:20
    %S(I)=sym('1');
   S(I)=V(I)*[(A(I+1)-X)^3]/(A(I+1)-A(I))/6+V(I+1)*[(X-A(I))^3]/(A(I+1)-A(I))/6+[B(I+1)-[V(I+1)*(A(I+1)-A(I))^2]/6]*(X-A(I))/(A(I+1)-A(I))+[B(I)-[V(I)*(A(I+1)-A(I))^2]/6]*(A(I+1)-X)/(A(I+1)-A(I));
end
%----------------------------------
Y='log(x)';
fplot(Y,[2,4],'r');%－－－擬合前函數的圖像
hold on
%-------------------樣條擬合曲線
for I=1:20
   ezplot(S(I),[A(I),A(I+1)]);
end
ezplot(F,[2,4]);%－－－擬合后函數的圖像  
title('函數ln(x)的Lagrangian以及三次樣條擬合圖像');
axis([2,4,0,2]);