##################################### # Find the inverse of a square matrix ##################################### # Define augmented matrix N:= 3: # Size of matrix to be inverted M:=2*N: A := array([[1,2,-1,1,0,0],[-4,1,1,0,1,0],[1,2,0,0,0,1]]); ### # Loop over columns. At each sweep, create unit pivot # and reduce entries underneath to zero by row operations ### for k from 1 to N do # Make 1st entry unit pivot div:= A[k,k]; for j from 1 to M do A[k,j]:= A[k,j]/div end do: # Display result print(A); # Perform row reductions for i from k+1 to N do # print("i=",i); q:= A[i,k]; for j from 1 to M do # print(j); A[i,j]:= A[i,j] - q*A[k,j] end do end do: # Display result print(A); end do: #### # Now matrix has 1's on leading diagonal. Transform to reduced echelon form #### print("Making into reduced echelon form..."); for k2 from 1 to N-1 do k:= N-k2+1; # reverse order # Perform row reductions for i2 from 1 to N-k2 do i:= (N-k2)-i2+1; # reverse order # print("i=",i); q:= A[i,k]; for j from 1 to M do # print(j); A[i,j]:= A[i,j] - q*A[k,j] end do end do: # Display result print(A); end do: # Pick out inverse from augmented form: B:= matrix(3,3,(i, j) -> A[i,j+3]); ##### # Using inbuilt commands ##### # with(linalg): # A1 := array([[1,2,-1],[-4,1,1],[1,2,0]]); # inverse(A1);