Lower/Upper Triangular Matracies and LU Factorization

[playboy]

Homework Statement

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.

Homework Equations

PLEASE IGNORE THIS PART FOR NOW AND SKIP TO

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)

n=length(A)
U1=A
L1=eye(n);
..for p=1:n
...for q=p+1:n
...L(q,p)=U(q,p)/U(p,p)
...U(q,p:n)=U(q,p:n)-U(p,p:n)*L(q,p)l
...end
..end



Forward Substition:

function x=forsub(L,b)

[n,n]=size(L);
x=zeros(n,1);
..for i=1:n
... x(i)=(b(i)-L(i,1:i-1)*x(1:i-1)/L(i,i)
..end


Backwards Substition:

function x=baksub(U,b)

[n,n]=size(L);
x=zeros(n,1);
..for i=n:-1:1
... x(i)=(b(i)-L(i,i+1:n)*x(i+1:n)/U(i,i)
..end

The Attempt at a Solution

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

But apart from trying to write the code, I am 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
and
Backword substitutions solves for x in Ux=b

I do not know where to go from here

Could somebody help me out please?

 

[AlephZero]

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.

 

[Mike_Li]

Hi,

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.

Cheers,

Mike Li.
New Zealand.