Borrar filtros
Borrar filtros

How to invert a square matrix using Gaussian elimination

26 visualizaciones (últimos 30 días)
Conor Pierce
Conor Pierce el 9 de Abr. de 2024
Editada: Torsten el 10 de Abr. de 2024
Hi there,
I have attached below the code I have written which is supposed to invert a square n x n matrix by Gaussian elimination. This m-file is as part of a larger project where it is implemented into a Newton-Raphson iterative method. The formula for this is Xn = Xn - (Df^-1)*(f). When i use this Gaussian code to invert the matrix, the Newton-Raphson method does not converge so any help on solving the problems in my Gaussia code attached below would be brilliant, thank you!
A = [1 1;2 5];
A_inv = Gauss_inverse(A)
A_inv = 2x2
0 1 1 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
A*A_inv
ans = 2x2
1 1 5 2
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
function A_inv = Gauss_inverse(A)
% Find dimensions of input matrix
[p,q] = size(A);
% If input matrix is not square, stop function
if p ~= q
error('Input must be a square n x n matrix')
end
% Target identity matrix to be transformed into the output
% inverse matrix
A_inv = eye(p);
%The following code actually performs the matrix inversion by working
% on each element of the input
for j = 1 : p
% Pivot row swapping for improved stability
[~,max_row] = max(abs(A(j:p, j))); % Find the row index of maximum element in the column
max_row = max_row + j - 1; % Adjust index to global reference
if A(max_row,j) ~= 0 % Skip if already zero
% Swap rows in both A and A_inv
temp = A(j,:);
A(j,:) = A(max_row,:);
A(max_row,:) = temp;
temp = A_inv(j,:);
A_inv(j,:) = A_inv(max_row,:);
A_inv(max_row,:) = temp;
end
end
if A(j,j) ~= 0
% Scale the pivot row to have 1 on diagonal
t = 1/A(j,j);
A(j,:) = t * A(j,:);
A_inv(j,:) = t * A_inv(j,:);
end
% Eliminate non-zero elements below the pivot
for i = j+1 : p
if A(i,j) ~= 0
t = -A(i,j);
A(i,:) = A(i,:) + t * A(j,:);
A_inv(i,:) = A_inv(i,:) + t * A_inv(j,:);
end
end
end

Respuestas (2)

Torsten
Torsten el 9 de Abr. de 2024
Editada: Torsten el 10 de Abr. de 2024
The newton update is
x_(n+1) = x_n - Df^(-1)*f,
not
x_(n+1) = x_n + Df^(-1)*f.
But you should not explicitly compute the inverse, but solve the linear system
Df*dx = -f
and set
x_(n+1) = x_n + dx
And as you can see from the simple example from above: your code to find the inverse of a matrix does not work as expected.

John D'Errico
John D'Errico el 9 de Abr. de 2024
Editada: John D'Errico el 9 de Abr. de 2024
I assume you are forced to use your own code for gaussian elimination here. A bad idea if you are not. NEVER write such code yourself, when high quality code written by professionals in the art is available.
If you are forced to write your own Newton method, An excellent suggestion is to do as Torsten suggests, to instead solve a system of equations. Best of course is to use backslash.
Ok, if you are forced to use an inverse there, and one you have written yourself, then first, verify that your code works! Never just throw something together and not test the pieces!
Once you have verified that piece of your code, then you need to recognize that Newton's method does not always converge. There is no assurance of convergence for all such problems. Sorry, but there is not. So, IF you are seeing non-convergence, you would want to learn to use the debugger. Are the Newton steps going in the correct direction? Does your code do something that makes sense?

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by