Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

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

Asked by Mark on 24 Nov 2011

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

Mark

Products

No products are associated with this question.

2 Answers

Answer by Walter Roberson on 24 Nov 2011
Accepted answer

http://www.mathworks.com/help/techdoc/ref/lu.html for LU factorization.

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

0 Comments

Walter Roberson
Answer by Derek O'Connor 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

http://www.derekroconnor.net/NA/LE/LE-2006-6.pdf

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).

3 Comments

Mark 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 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 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.

Derek O'Connor

Contact us