求矩阵函数.m

计算任意矩阵函数的matlab程序

%此文件用于计算任意矩阵函数的值，在程序运行前，要先设定矩阵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    %最后的结果