%此文件用于計算任意矩陣函數的值，在程序運行前，要先設定矩陣A的內容和矩陣函數的形式
A=[4 6 0;-3 -5 0;-3 -6 1]; %矩陣函數中的矩陣
f=sym('sin(x)' );  %用符號標識函數形式

%計算矩陣特征值
lamda=(eig(A))';%lamda表示特征值
n=length(lamda);
lamda_value=[];
lamda_jieshu=[];

%求出相異特征值和其重數
while length(lamda)>0
    identity=find(abs(lamda-lamda(1))<0.01);  %尋找相同的特征值
    lamda_value=[lamda_value,lamda(1)];
    lamda_jieshu=[lamda_jieshu,length(identity)];  %階數指的是相同的特征值的個數
    lamda(identity)=[]; %去掉相同的特征值項
end
lamda_count=length(lamda_value);   %不同的特征值的個數

%用牛頓插值求r(z)
newton=zeros(n,n);
startline=0;     %用于標志行的起始位置
for k=1:lamda_count
    if k==1
        startline=1;
    else
        startline=startline+lamda_jieshu(k-1);  % 用于不同特征值的處理，k代表每個特征值
    end
    
    for l=1:lamda_jieshu(k)   %開始處理每個特征值，l表示相同特征值中的第幾個
        currentline=startline+(l-1);  %currentline指的是當前處理的行數，即當前在處理第幾個特征值
        newton(currentline,1)=subs(f,lamda_value(k));  %newton矩陣是要求的差值矩陣，subs函數的作用是替換，即講newton中的元素用lamda_value(k)經f處理后的值替代
        if currentline>=3      %處理newton矩陣中下三角中對角線以下的元素
            for m=2:currentline-1
                if l==1
                    newton(currentline,m)=(newton(currentline,m-1)-newton(currentline-1,m-1))/(lamda_value(k)-lamda_value(k-1));%如果是第一個，則按照正常的插值公式計算
                else
                    newton(currentline,m)=subs(diff(f,m-1),lamda_value(k))/factorial(m-1); %由于是相同的特征值，則求出的結果是導數的形式，其中factorial函數是取階乘的作用
                end
            end
        end
        if currentline>=2   %處理newton矩陣中對角線上的元素
            if k~=1  %如果是第一個具有多個相同特征值的元素，則進行求導處理
                newton(currentline,currentline)=(newton(currentline,currentline-1)-newton(currentline-1,currentline-1))/(lamda_value(k)-lamda_value(1));%直接處理
            else
                newton(currentline,currentline)=(subs(diff(f,currentline-1),lamda_value(k)))/factorial(currentline-1); %若不是第一行，則由于是相同的元素，故結果是導數的形式,diff函數用于求導
            end
        end
    end
end

%計算得到多項式的表達式
syms z;     %定義一個變量z
nf=0;       %用于計算多項式的形式
nf_part=1;  %多項式的一部分
for k=1:lamda_count   
    if k==1
        startline=1;
    else
        startline=startline+lamda_jieshu(k-1);
    end
    for l=1:lamda_jieshu(k)
        currentline=startline+(l-1);
        if l==1
            if k~=1
                nf_part=nf_part*(z-lamda_value(k-1)); %乘以上一個特征值構成的多項式因子
            end
        else
            nf_part=nf_part*(z-lamda_value(k)); %重復乘以特征值構成的多項式因子
        end
        if currentline>1
            nf=nf+nf_part*newton(currentline,currentline); % 乘以系數，即乘以插值
        end
    end
end
nf=nf+newton(1,1); %得到最后的多項式形式

%將A代入求出nf(A),即得該矩陣函數的值
nf_coff=sym2poly(nf);  %獲取多項式系數，sym2poly函數用于返回多項式的系數
result=zeros(n,n);
for k=length(nf_coff):-1:1
    result=result+nf_coff(k)*A^(length(nf_coff)-k); %計算最后的結果
end
result    %最后的結果