Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Lower/Upper Triangular Matracies and LU Factorization

  1. Feb 4, 2007 #1
    1. The problem statement, all variables and given/known data

    Let A and B be invertible n×n matrices and b be an n×1 vector.

    Write a MATLAB function with inputs (A,B, b) to solve the equation x=B^−1*(2A^−1 + 1)b

    Make use of functions "LU Facotrization, Forward Substituion and Backwards Substitution" and DO NOT calculat any matrix inverse.

    2. Relevant equations

    PLEASE IGNORE THIS PART FOR NOW AND SKIP TO 3. The attempt at a solution

    Well, I came up with the functions for LU Facotrization, Forward Substituion and Backwards Substitution;

    LU Facotrization:

    function [L,U]=lufactor(A)

    ..for p=1:n
    ....for q=p+1:n

    Forward Substition:

    function x=forsub(L,b)

    ..for i=1:n
    ..... x(i)=(b(i)-L(i,1:i-1)*x(1:i-1)/L(i,i)

    Backwards Substition:

    function x=baksub(U,b)

    ..for i=n:-1:1
    ..... x(i)=(b(i)-L(i,i+1:n)*x(i+1:n)/U(i,i)

    3. The attempt at a solution

    I wrote the above functions, so that is a start,

    But apart from trying to write the code, Im kind of lost in understanding the math - which is much more important.

    I want to LU-factorize A and B to get and upper and lower matrix for each of A and B.

    With these matracies, how would I go abouts get it's inverse?

    Forword substitution solves for x in Lx = b
    Backword substitutions solves for x in Ux=b

    I do not know where to go from here

    Could somebody help me out please?
  2. jcsd
  3. Feb 5, 2007 #2


    User Avatar
    Science Advisor
    Homework Helper

    To solve Ax = b, first find L and U. Then

    Ax = b
    (LU)x = b

    Since matrix multiplication is associative you can change the position of the brackets:

    L(Ux) = b

    Ux is a vector - call it y

    First you solve

    Ly = b (forward substitution)

    Then solve

    Ux = y (backward substitution).

    BTW I didn't check your code in detail, but you are doing the right sort of things.

    To make up an "easy" test problem, you can work backwards: write down some "simple" matrices L and U (say 3x3), then multiply them to get A, then write down a "simple" answer x, then multiply Ax to get b.

    Then, test if your functions can work back from A and b to find the L U and x that you first thought of.
  4. Jun 18, 2009 #3

    A really neat way - and outside of MatLab a really fast way - of solving the inverse of a matrix is to recognise that any matrix, times it's inverse, = the identity matrix.

    Which means that when you have a decomposition of any variety - L - then you find the inverse of L to begin with using only forward substitution as follows:

    L.Li = I

    Where I is the identity matrix and Li is the inverse of L

    This is solved by solving 1 column of the inverse at a time:

    L.Li = [1,0,0,0.....]

    and storing that column, then looping back and solving with [0,1,0,0,0.....0], then [0,0,1,0,0,0....0] & etc.

    The full inverse can then be computed as Li.Li' - i.e. the inverse of the triangular decomposition (or factor) multiplied by it's own transpose.

    If you use CSR to store the data, it's really fast too - except for the Triangular x transpose at the end, which can't be avoided.


    Mike Li.
    New Zealand.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook