MATLAB Answers

## Solve for x in (A^k)*x=b (sequentially, LU factorization)

Asked by Mark

### Mark (view profile)

on 24 Nov 2011
Latest activity Edited by Rena Berman

### Rena Berman (view profile)

on 9 Jul 2018
Accepted Answer by Walter Roberson

### Walter Roberson (view profile)

Without computing A^k, solve for x in (A^k)*x=b.
A) Sequentially? (Pseudocode)
for n=1:k
x=A\b;
b=x;
end
Is the above process correct?
B) LU factorizaion?
How is this accompished?

#### 0 Comments

Sign in to comment.

## 2 Answers ### Walter Roberson (view profile)

Answer by Walter Roberson

### Walter Roberson (view profile)

on 24 Nov 2011
Accepted Answer

However, I would suggest that LU will not help much. See instead http://www.maths.lse.ac.uk/Personal/martin/fme4a.pdf

Nicholas Lamm

### Nicholas Lamm (view profile)

on 9 Jul 2018
A) Linking to the documentation is about the least helpful thing you can do and B) youre not even right, LU decomposition is great for solving matrices and is even cheaper in certain situations.

Sign in to comment.

Answer by Derek O'Connor

### Derek O'Connor (view profile)

on 28 Nov 2011

Contrary to what Walter says, LU Decomposition is a great help in this problem. See my solution notes to Lab Exercise 6 --- LU Decomposition and Matrix Powers
Additional Information
Here is the Golub-Van Loan Algorithm for solving (A^k)*x = b
[L,U,P] = lu(A);
for m = 1:k
y = L\(P*b);
x = U\y;
b = x;
end
Matlab's backslash operator "\" is clever enough to figure out that y = L\(P*b) is forward substitution, while x = U\y is back substitution, each of which requires O(n^2) work.
Total amount of work is: O(n^3) + k*O(n^2) = O(n^3 + k*n^2)
If k << n then this total is effectively O(n^3).

Mark

### Mark (view profile)

on 28 Nov 2011
Thank you very much. Still, I'm not quite sure how to solve for y in Ly=Pb, y=inv(L)*P*b? If so, what about x=invU)*(inv(L)*P*b)? Last, b=inv(U)*(inv(L)*P*b) and then return b in place of x?
Derek O'Connor

### Derek O'Connor (view profile)

on 28 Nov 2011
Once you have L, U, and P, then Ax = b is replaced by LUx = Pb. This equation is solved in 2 steps, where y = Ux and c = Pb:
1. Solve Ly = c for y using forward substitution, because L is lower triangular.
2. Solve Ux = y for x using back substitution, because U is upper triangular.
Note y = L^(-1)c and x = U^(-1)y are mathematical statements and are never used to numerically solve these equations.
I think you may need to do some reading and thinking about LU factorization(decomposition).
Take a look at my webpage http://www.derekroconnor.net/NA/na2col.html for a list of textbooks and Section 5 Notes on Numerical Linear Algebra, pages 12, 13, and 21-25.
Or, better still, read the relevant chapter in Cleve Moler's book which is available online at http://www.mathworks.com/moler/
Derek O'Connor

### Derek O'Connor (view profile)

on 28 Nov 2011
Oh dear. It has just struck me that this may be a homework problem and I have given the game away.

Sign in to comment.