Precise Matrix Power: A^k

Hello,

I am currently fighting with (more than) double-precision numerics, with which I am not familiar, so the answer probably lies within this topic. I would like to calculate

with and but obviously the precision of will be bad for big k. Is there a datatype of which I am not aware, or any other trick that will help me? Thanks in advance. :)

Answer 1

[Rik]

It will probably be helpfull to use the Symbolic Math Toolbox, and use the tips explained here. For mod(a^b,k) there is a direct function, but for your use case there doesn't seem to exist one.

[Eric Schöneberg]

I've generated the following function:

function A = matrixPower(Matrix, exponential, precision)
A = vpa(eye(size(Matrix)), precision);
for k=1:exponential
   A = vpa(A*Matrix, precision); 
end
A = double(A);
end

which seems to work. I've tested the method using this script:

A = magic(5); 
precision = 100;
for k = 0:50
   B = matrixPower(A, k, precision);
   C = A^k;
   Error = B-C
end

and the Error(-Matrix) was Zero for k=1:10 and grew to 1.0e+74*Matrix for k = 50. Thanks again for your input @Rik.