Why does the computational error of this matrix exponential algorithm has this shape?

113 Views Asked by At

enter image description here

Matrix exponential is defined as: enter image description here

Why is that the computational error decreases until around k = 75, then it stays the same? Is it because the factorial is too big, and because Matlab only uses 16 digits of precision, the round-off error is neglected?

Here's my Matlab algorithm:

format long
%     load('CA3matrix.mat');
    A = randn(500);
    terms_arr = 5:5:150; % number of terms in the exp(A)
    prev_term = 1; % store the previous terms, avoid duplicate calculation
    expAk = A^0;
    errors = [];

    for i = terms_arr
       for k = 1:i
            term = (1/k) * A * prev_term;
            prev_term = term;
            expAk = expAk + term;
            %expAk = expAk + (1/factorial(k)) * A^(k);
       end
       expA = expm(A);
       err = norm(expA - expAk)/norm(expA);
       errors = [errors err];
       expAk = A^0;
       prev_term = 1;
    end

    semilogy(terms_arr, errors, '*')
    title('Computational error versus k')
    xlabel('k','fontsize',16)
    ylabel('Computational error','fontsize',16)
    set(gca,'FontSize',14)