Inverse Matrix script is breaking with anything more than 3x3 matrix.

Hi,

For uni we need to build a script that can find the inverse of any nxn matrix, using eigenvalues to do this.

clear
%User sets up matrix with nxn dimensions. 
A=input('Input matrix A: ');
%Aquire EigenVlaues/Vectors
[V,D]=eig(A);
%Process eigs to create inverse of A
Ainv = V*(inv(D))*V.'

This is my code so far, which works for 3x3 and below. It uses the equation av=λv to find the inverse of a.

However 4x4 and above it breaks. Can anyone explain why it breaks and also how to make it work for any nxn matrix?

Thanks in advance.

[John D'Errico]

First of all, you use the function inv inside your algorithm! If that was acceptable, then a perfectly good solution would be as simple as...

clear
%User sets up matrix with nxn dimensions. 
A=input('Input matrix A: ');
%Aquire EigenVlaues/Vectors
[V,D]=eig(A);
%Process eigs to create inverse of A
Ainv = inv(A) + 0*V + 0*D;

You must agree that this will return the inverse of A whenever that inverse exists, AND that it uses both of the outputs of eig?

Surely there is a fundamental problem with that idea. You are supposed to use eig to solve for the inverse of a matrix, in a way that does not just use inv. But if you then use inv inside, then it seems a bit silly. And it would probably earn you a poor grade on your assignment.

So you might consider if there is a way to find the inverse of a diagonal matrix, yet does NOT use inv? (Is D always diagonal?) How would you form the inverse of such a matrix? What does the function diag do?

A for why it fails, what evidence do you have that it truly does fail? Why do you think that it fails, yet you think it works for (any?) 3x3 matrix? Just for kicks...

A = rand(3);
[V,D] = eig(A);
inv(A) - V*inv(D)*V.'
ans =
       37.356      -54.806      -4.9635
       36.887      -52.593      -8.2745
      -85.759       124.12       15.237

Hmm. If it really did work on any 3x3 matrix, then why is that not all zeros, or at least close to zero?

Even better, try this repeatable example:

A = [1 1 1;0 1 1;0 0 1]
A =
     1     1     1
     0     1     1
     0     0     1

[V,D]= eig(A)
V =
            1           -1            1
            0   2.2204e-16  -2.2204e-16
            0            0   4.9304e-32
D =
     1     0     0
     0     1     0
     0     0     1

inv(A)
ans =
     1    -1     0
     0     1    -1
     0     0     1

V*inv(D)*V.'
ans =
            3  -4.4409e-16   4.9304e-32
  -4.4409e-16   9.8608e-32  -1.0948e-47
   4.9304e-32  -1.0948e-47   2.4309e-63

Hmm. That would seem to be a problem.

So one idea might be to find a solution that will work at least on symmetric non-singular matrices, that does not directly use inv.