1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: Maple - LU Factorization with Partial Pivoting

  1. Jan 12, 2012 #1
    1. Im asked to write a procedure that will inverse square matices using LU factorization with partial pivoting.

    2. Im also told that the procedure should return the inverse matrix and report an error if it cannot do so.

    3. So far i've come up with the code below but am not confident it is correct and I am also really struggling to find help with the next step for solving the linear systems i need too proceed with. Any adivce or suggestions:
    restart :

    n:= n; # all matrices are nxn matrices
    A:= array (1..n, 1..n);
    L:= array (1..n, 1..n);
    U:= array (1..n, 1..n);
    LU:= array (1..n, 1..n);

    print('A is', A);

    # carry out the Gaussian elimination
    for j from 1 to n-1 do
    for i from j+1 to n do
    # calculate the multiplier for the ith row at the jth stage
    m:= A[i,j]/A[j,j];
    A[i,j]:= m;
    for k from j+1 to n do
    A[i,k]:= A[i,k]-m*A[j,k];
    print ('new A is', A);

    # Extract L and U from the new A
    for i from 1 to n do
    for k from 1 to n do
    if i<k then
    U[i,k]:= A[i,k]
    L[i,k]:= 0;
    elif i=k then
    U[i,k]:= A[i,j]
    U[i,k]:= 0
    L[i,k]:= A[i,j]

    print('L is',L);
    print('U is',U);

    # verify that A=LU
    for i from 1 to n do
    for k from 1 to n do
    LU[i,k]:= add(L[i,j]*U[i,k], k=1..n);

    print('LU is', LU);
    1. The problem statement, all variables and given/known data

    2. Relevant equations

    3. The attempt at a solution
  2. jcsd
  3. Jan 12, 2012 #2
    So, i've been working on this problem for a few hours today and have come up with this procedure (below) but need advice on how to change it so that it works for a square matrix only and then finally how to form the inverse matrix to finish the question. Any ideas welcome...
    # LU decomposition using Gaussian elimination with partial pivoting
    # lufactor does the LU factorization
    # lusolve solves a linear system using the LU factorization
    lufactor := proc(A)
    local n, r, k, p, i, temp, m, j;

    n := rowdim(A); # get number of rows in matrix A
    r := vector(n); # create pivot vector
    for k from 1 to n do r[k] := k od; # initialize pivot vector

    for k from 1 to n-1 do # loop over columns

    # We find the maximal element in absolute value among the
    # entries A[r[k],k], A[r[k+1],k], .... A[r[n],k].
    # Interchange entries in pivot vector

    p := k;
    for i from k+1 to n do
    if (abs(A[r,k]) > abs(A[r[p],k])) then p := i; fi;
    if A[r[p],k]=0 then
    ERROR(`Matrix is singular`); # procedure will exit here
    if (p <> k) then # interchange entries k and p in pivot vector
    temp := r[k];
    r[k] := r[p];
    r[p] := temp;

    # now do a gaussian elimination step

    for i from k+1 to n do # loop down column below diagonal
    m := A[r,k]/A[r[k],k]; # get multiplier
    for j from k+1 to n do # subtract multiple of row i from
    A[r,j] := A[r,j] - m*A[r[k],j]; # row k
    A[r,k] := m; # save multiplier
    RETURN(op(r)); # return pivot vector

    # solve the linear system Ax=b given the LU factorization in a returned
    # from ludecomp, the pivot vector r returned from ludecomp and the right
    # hand side vector b.

    lusolve := proc(A,r,b)
    local n, y, x, i, s, j;
    n := rowdim(A); # get number of rows in matrix A
    y := vector(n); # create vector for solution of Ly=Pb
    x := vector(n); # create vector for solution of Ux=y

    # do forward substitution to solve Ly=Pb

    y[1] := b[r[1]];
    for i from 2 to n do
    s := b[r];
    for j from 1 to i-1 do
    s := s - A[r,j]*y[j];
    y := s;

    # do backward substitution to solve Ux=y

    x[n] := y[n]/A[r[n],n];
    for i from n-1 to 1 by -1 do
    s := y;
    for j from i+1 to n do
    s := s - A[r,j]*x[j];
    x := s / A[r,i];
    RETURN(op(x)); # return solution vector
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook