- #1
Rupert11
- 2
- 0
1. I am asked to write a procedure that will inverse square matices using LU factorization with partial pivoting.
2. I am 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];
od:
od:
od:
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]
L[i,k]:=1;
else
U[i,k]:= 0
L[i,k]:= A[i,j]
fi;
od:
od:
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);
od:
od:
print('LU is', LU);
2. I am 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];
od:
od:
od:
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]
L[i,k]:=1;
else
U[i,k]:= 0
L[i,k]:= A[i,j]
fi;
od:
od:
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);
od:
od:
print('LU is', LU);