matrixfun.m

来自「Matlab程序」· M 代码 · 共 26 行

function MatrixB=matrixfun(Coef,MatrixA)
syms x;
f=feval(Coef,x);
EigA=sort(eig(MatrixA));
LengthEigA=length(EigA);
NewtenMatrix=zeros(LengthEigA);
NewtenMatrix(:,1)=subs(f,x,EigA);

for i=2:LengthEigA
    for j=i:LengthEigA
        if abs(EigA(j)-EigA(j-1))<1e-14;
            NewtenMatrix(j,i)=subs(diff(f,i-1,x),x,EigA(j))/prod(linspace(1,i-1,i-1));
        else
            NewtenMatrix(j,i)=(NewtenMatrix(j,i-1)-NewtenMatrix(j-1,i-1))/(EigA(j)-EigA(j-i+1));
        end
    end
end

MatrixB=0;
for i=LengthEigA:-1:1
    if i==LengthEigA
       MatrixB=NewtenMatrix(i,i)*eye(LengthEigA);
   else
       MatrixB=NewtenMatrix(i,i)*eye(LengthEigA)+(MatrixA-EigA(i)*eye(LengthEigA))*MatrixB;
   end
end