how can I calculate A^n when n is a symbolic positive integer?

A = [1/2,1/2,0,0;1/2,0,1/2,0;1/2,0,0,1/2;0,0,0,1]
syms n positive integer
A^n

This would not work because it got stuck and could never stop.

Answer 1

[Walter Roberson]

syms N positive integer 
[V, D] = eig(sym(A)) 
result = V*diag(diag(D).^N)/V

Note those are the matrix operations * and / not element by element operations

[Songbai Jin]

Thank you for your help.

I find

result = V*diag(diag(D).^N)*inv(V)

works a little better than '/' operation.

However, it seems MATLAB is not good at doing symbolic calculation. The result is too complex and always contains imaginary unit i, which should not appear in the result.

[Songbai Jin]

This method does not work if A is NOT similar to diagonal matrix.

[Walter Roberson]

A = [1/2,1/2,0,0;1/2,0,1/2,0;1/2,0,0,1/2;0,0,0,1];
syms N positive integer 
[V, D] = eig(sym(A)); 
result = V*diag(diag(D).^N)/V;
result5 = subs(result,N,5);
direct_version = sym(A)^5;

vpa(real(result5),10) - vpa(direct_version,10)

ans = 

[Walter Roberson]

The method does require that the matrix is diagonizable.

A = {[1/2,1/2,0,0;1/2,0,1/2,0;1/2,0,0,1/2;0,0,0,1], [1 3; 0 1]}

A = 1×2 cell array
    {4×4 double}    {2×2 double}

format long g
syms N positive integer 
for K = 1 : numel(A)
    [V, D] = eig(sym(A{K})); 
    if rank(V) ~= size(A{K},1)
        fprintf('diagonalization failed for matrix #%d\n', K);
    else
        result = V*diag(diag(D).^N)/V;
        result5 = subs(result,N,5);
        direct_version = sym(A{K})^5;
        fprintf('matrix #%d\n', K)
        double(result5 - direct_version)
    end
end

matrix #1

ans = 
       5.94278939548003e-71 -  6.91910479616603e-71i      -3.22608567183202e-71 +  1.17087100589517e-71i      -8.34820415079338e-72 -  6.46552480659025e-72i        7.3577392515467e-72 +  6.80939948034781e-75i
       2.16204338006988e-70 -  8.95310060594234e-71i       1.98092979849334e-71 +  1.79094167932363e-70i      -7.69732721700271e-71 -  2.37004100891168e-72i                          0 +  5.48748524668485e-75i
      -5.34851045593203e-71 -  2.10120053625901e-71i      -2.23562077258534e-71 +  4.81259819794677e-71i       5.23531446744669e-72 +  1.01119582283165e-72i       1.58474383879467e-71 -  5.10704961026086e-75i
                          0 +                     0i                          0 +                     0i                          0 +                     0i                          0 +                     0i

diagonalization failed for matrix #2

At the moment, I am not aware of any general formula for cases that are not diagonalizable.

[Walter Roberson]

In cases where the matrix is numeric, you can approximate by using numeric eigenvalues, which removes problems with working with matrices 5 x 5 or larger (where it might not be possible to algebraically find the exact eigenvalues)

A = {[1/2,1/2,0,0;1/2,0,1/2,0;1/2,0,0,1/2;0,0,0,1], [1 3; 0 1]}

A = 1×2 cell array
    {4×4 double}    {2×2 double}

format long g
syms N positive integer 
for K = 1 : numel(A)
    [V, D] = eig((A{K})); 
    if rank(V) ~= size(A{K},1)
        fprintf('diagonalization failed for matrix #%d\n', K);
    else
        sV = sym(V);
        result = sV*diag(diag(D).^N)/sV;
        result5 = subs(result,N,5);
        direct_version = (A{K})^5;
        fprintf('matrix #%d\n', K)
        double(result5 - direct_version)
    end
end

matrix #1

ans = 
      -1.63526849241225e-16 - 2.96770310333167e-167i       -1.4643637431556e-17 + 2.47308591944306e-167i      -1.81451191767119e-16 - 3.09769808195142e-168i       3.59621678439901e-16 +                     0i
      -1.83427542145719e-16 + 5.72773437856958e-168i      -4.98395007893118e-17 + 2.38476142232009e-167i      -1.29593823909905e-16 -  1.0598939654756e-168i       3.62860866844936e-16 +                     0i
       -1.0923608999026e-16 - 1.69583034476095e-167i      -4.03961279131375e-17 - 1.48385155166583e-167i      -9.40591452819007e-17 +                     0i       2.43691363185299e-16 + 1.41319195396746e-168i
                          0 +                     0i                          0 +                     0i                          0 +                     0i                          0 +                     0i

diagonalization failed for matrix #2

[Walter Roberson]

The result is too complex and always contains imaginary unit i, which should not appear in the result.

Let us have a look at determining the roots of a real-valued polynomial of degree 3 that has real roots

C = [1 2 -3 -4]

C = 1×4

1 2 -3 -4

syms a3 a2 a1 a0 x
p = a3*x^3 + a2*x^2 + a1*x + a0

p = 

sols = solve(p, 'maxdegree', 3)

sols = 

The above is the general solution, including the case where the coefficients might be complex valued and where the roots might be complex valued. Let us narrow it down by substituting in real coefficients:

subsols = subs(sols, [a3 a2 a1 a0], C)

subsols = 

Okay, so you would not necessarily expect the imaginary components to vanish immediately upon substituting in real coeficients. Let us simplify() the expressions; that should the trick:

ssubsols = simplify(subsols)

ssubsols = 

No, not good enough.

Ummm, maybe we need to expand the expressions and simplify the results of that?

sss = simplify(expand(ssubsols))

sss = 

Nope, there are still imaginary components in that. We have a cube root of a complex-valued expression. Do we have an reason to expect to be able to rewrite that (19+6*sqrt(51)*i)^(1/3) in a separable form A+B*i so that we would have a hope of being able to rationalize the expressions to eventually symbolically get rid of all of the *i components? I do not think so.

What are the roots?

double(sss)

ans = 3×1

1.5616 -1.0000 -2.5616

Or in the more direct form,

solve(poly2sym(C))

ans = 

double(ans)

ans = 3×1

-1.0000 -2.5616 1.5616

So the roots are real-valued... but can we really expect to create a generalized formula that does not involve i anywhere that would return the real roots just given the coefficients and the promise that the roots are all real?

[Songbai Jin]

Thank you for your help! It really helps me a lot!

Although I do not know how to remove i from A^N result, I find a way to get real root result.

C = [1 2 -3 -4];
syms a3 a2 a1 a0 x
p = a3*x^3 + a2*x^2 + a1*x + a0;
sols = solve(p, 'maxdegree', 3);
subsols = subs(sols, [a3 a2 a1 a0], C)

subsols = 

result = simplify(subsols,'Steps',100)

result = 

This method does not work in A^N result, probably because it needs more steps, or MATLAB just can not simplify the A^N answer.

[Songbai Jin]

What's more, I find the mpower function implemented by MATLAB needs improving.

A = [1,3;0,1];
syms N positive integer
A^N

ans = 

A = [1,1;0,0];
syms N positive integer
A^N

Error using ^
Singularity.

The function throws error if A is singular matrix, but it can calculate A^N if A is not diagonalizable.

Answer 2

[KSSV]

A = [1/2,1/2,0,0;1/2,0,1/2,0;1/2,0,0,1/2;0,0,0,1]

A = 4×4

0.5000 0.5000 0 0 0.5000 0 0.5000 0 0.5000 0 0 0.5000 0 0 0 1.0000

syms n positive integer
A.^n

ans = 

[Dyuman Joshi]

That wouldn't be necessarily same as A^n

[Songbai Jin]

I hope to calculate matrix multiplication A^n not elementwise product A.^n.

Thank you for your help.

[KSSV]

A = [1/2,1/2,0,0;1/2,0,1/2,0;1/2,0,0,1/2;0,0,0,1]
n = 10 ; 
B = eye(size(A)) ; 
for i = 1:n
    B = B*A ; 
end