Trying to write a program that calculates the inverse of a 3x3 matrix. My program works for some matrices, but not for all. Can someone please look at my code and assist me?

%calculate determinant 
m = input('Enter a 3x3 matrix using bracket notation: ');
a = m(1,1);
b = m(1,2);
c = m(1,3);
d = m(2,1);
e = m(2,2);
f = m(2,3);
g = m(3,1);
h = m(3,2);
i = m(3,3);
detA = a*(e*i-f*h);
detB = b*(d*i-f*g);
detC = c*(d*h-e*g);
det = detA - detB + detC;

%calculate cofactor
a2 = e*i-h*f;
b2 = -d*i-g*f;
c2 = d*h-g*e;
d2 = -(b*i-h*c);
e2 = a*i-g*c;
f2 = -(a*h-g*b);
g2 = b*f-e*c; 
h2 = -(a*f-d*c);
i2 = a*e-d*b;
detA2 = a2*(e2*i2-f2*h2);
detB2 = b2*(d2*i2-f2*g2);
detC2 = c2*(d2*h2-e2*g2);
cof = detA2 - detB2 + detC2;
%cofactor matrix
n = [a2 b2 c2; d2 e2 f2; g2 h2 i2];

%calculate transpose of cofactor
v(:,1) = n(1,:);
v(:,2) = n(2,:);
v(:,3) = n(3,:);

%calculate inverse
if(det ~= 0)
    inv = (1/det)*v;
    disp(inv);
else
    fprintf('Matrix is not invertible.\n');
end;

[Jacky Jo]

In matlab we don't have to do as you did. Try the following:

c

lc; clear all; close all;
m=input('Enter a 3x3 matrix using bracket notation: \n'); 
Det=det(m)

if(Det ~= 0)
    format bank;
    fprintf('\nThe determinant is:\n  %d\t%d\t%d\n %d\t%d\t%d\n %d\t%d\t%d\n',inv(m));  
else
    fprintf('\nMatrix is not invertible!\n');
end;

[William Diuguid]

I can't use any MATLAB functions. I have to do it myself

[Walter Roberson]

http://uk.mathworks.com/matlabcentral/answers/38787-what-can-be-programmed-without-any-built-in-functions

[Walter Roberson]

Give us an example of a matrix it fails on.

[William Diuguid]

[1 3 5; 5 3 1; -5 -2 3]. I checked online and the inverses aren't the same

[Walter Roberson]

When I run you code, the answer is the same as MATLAB's inv() routine to within round-off error.

By the way, I really advise that you find a variable name other than "inv" to use, as that is the name of the MATLAB inv() routine. Anything you can do to save wear and tear on the minds of people reading the code is good.

[William Diuguid]

I just ran my program with [1 9 8; 3 2 1; 4 5 6] and inv() with the same matrix and got different results

[Walter Roberson]

I get exactly the same values.

Enter a 3x3 matrix using bracket notation:  [1 9 8; 3 2 1; 4 5 6]
        -0.111111111111111         0.222222222222222         0.111111111111111
         0.349206349206349         0.412698412698413        -0.365079365079365
        -0.111111111111111        -0.492063492063492         0.396825396825397

>> clear inv   %badly named variable!
>> inv(m)

ans =

          -0.111111111111111         0.222222222222222         0.111111111111111
           0.222222222222222         0.412698412698413        -0.365079365079365
          -0.111111111111111        -0.492063492063492         0.396825396825397

[dpb]

Walter, ans(2,1) isn't within roundoff altho rest are...a quick glance might overlook that.

Answer 1

[Zeus]

This code works for any square and invertible matrix of any order

function A_inverse=inverse(A)
% A_inverse as the name suggests is the inverse of A
% A is a square matrix which is invertible
% if A is not a square matrix or A is square matrix but not invertible,then the output is not equal to inverse of A
a=length(A); % the order of A is axa
I=eye(a);
augmat=[A I]; % AUGMENTED MATRIX
% GAUSSIAN ELMINATION METHOD:

% when A is invertible, [A I] is row equivalent to [I inv(A)]
% in other words, the row operations that convert A to I also converts I to inv(A)
% I is identity matrix of order axa, inv(A) is the inverse of A of order axa

% Converting A to its Echelon form
for i=1:a-1
    m=augmat(i,i);
    augmat(i,:)=augmat(i,:)/m; % normalization,so that pivot values will be equal to 1
    for j=i:a-1
        augmat(j+1,:)=augmat(j+1,:) - augmat(i,:)*augmat(j+1,i); % making the elements below the pivots as zero
    end
end
augmat(a,:)=augmat(a,:)/augmat(a,a); % noralization of the last row of A

% Converting A from its Echelon form to Row Reduced Echelon form
for k=2:a
    for g=(k-1):-1:1
        augmat(g,:)=augmat(g,:)-augmat(k,:)*augmat(g,k); % makes the elements above pivots to be row
    end
end
%We end up with A converted to I and I will be converted to inv(A)
A_inverse=augmat(:,a+1:2*a); % extracting inv(A) from augmented matrix [I inv(A)]
end

Answer 2

[the cyclist]

Your definition of b2 has a sign error.