Thank you for all of your comments, John and Jan. Afterall I have decided to continue instead of giving up. Here, I want to say thank again to the encouragement from Jan.
The presentation in File Exchange is now updated to include two files: inv0.m and inv0det_pvt.m.
(1) AI = inv0(AO) -- Find AI(inverse) directly from AO(given matrix) without any pivoting. The code is only 4 statement lines total.
(2) [AI,det] = inv0det_pvt(AO,pvt) -- Find both AI and det(determinant) from AO, with or without pivoting. The "write-in" is used for the pivoting (pvt = 2). It is only good for AO of small order and not practical at all. However it may resolve the problem when the divisor 'd' in the k-loop is zero or close to zero. Such as d=AO(1,1)=0.
The option of automatic pivoting with the preset critique is the way to proceed. I think the critique for pivoting is to select the pivoting element such that the divisor d is the absolute maximum among all computed d's. If the value of d, resulted from the selected pivoting element, is zero or close to zero, then the given matrix AO must be singular.
Anyway the derived algorithm for matrix inversion is very pratical compare to some other classical approaches, such as Gauss-Jordan Elimination. I hope our FEX users will like to use it.
I am not majored in mathematics. I am retired now and soon approach to 80. I need some help regarding the option for automatic pivoting.
Thank you again for your comments. Afterward I do not expect people would want to use my code at all, since the MATLAB functions 'inv' and 'det' have already been here to be used conveniently. However, it does not show how these built-in functions are developed, such as by co-factors, eigenvalues, Guanssian elimination, etc. Perpaps, this is why people like to develop their own codes.
Sometimes, people want to derive his own routine that may be better than the built-in one in some extent. For example, the routine 'polyroots' presented in Mathworks Files Exchange is
more suitable than the built-in function 'roots' for any given polynomials with many high multiple roots.
If anyone like to develop his own codes about the inverse and determinant of a given square matrix, he may do it by applying the approach used in my code. It is very siple and straightforward, and only six statement lines.
Finally, the errors have been found during repeatly running the code
with pvt = 1 option. It will yield the determinant of correct
absolute value but with different sign, either + and - . It needs to be corrected.
I may withdraw this contribution if I feel no good.
I do appreciate John's comments, so that I am to force myself to do some further homework. Hope it will give him some answers but not expect all.
The code is now updated to include the pivoting scheme. There are three options to choose:
(1) pvt=0 ---> no pivoting,
(2) pvt=1 ---> pivoting elements randomly selected,
(3) pvt=2 ---> pivoting elements by writing-in.
If the very first element of the given square matrix is zero, it does surely fail for option(1), but it will be OK by sucessively running either option(2)or(3) for any non-singular matrix.
The code derived is very short (10 lines for the original and less than 30 for the updated). Yet it will give both inverse and determinant of a given square matrix.
The correction works for the reported problem, but, there is more. For instance, for p=[3e-04 1e-09 1e-09 5e-08 2e+04] results are wrong and in addition multiplicities appear as complex numbers.
There seems to be a bug when it comes to linear polynomials (which of course can not have multiple zeros). poly_roots([1 2]) does not return a zero and poly_roots([1 1 0 0 0]) returns just the triple zero 0.
Thank you for all of your comments, John and Jan. Afterall I have decided to continue instead of giving up. Here, I want to say thank again to the encouragement from Jan.
The presentation in File Exchange is now updated to include two files: inv0.m and inv0det_pvt.m.
(1) AI = inv0(AO) -- Find AI(inverse) directly from AO(given matrix) without any pivoting. The code is only 4 statement lines total.
(2) [AI,det] = inv0det_pvt(AO,pvt) -- Find both AI and det(determinant) from AO, with or without pivoting. The "write-in" is used for the pivoting (pvt = 2). It is only good for AO of small order and not practical at all. However it may resolve the problem when the divisor 'd' in the k-loop is zero or close to zero. Such as d=AO(1,1)=0.
The option of automatic pivoting with the preset critique is the way to proceed. I think the critique for pivoting is to select the pivoting element such that the divisor d is the absolute maximum among all computed d's. If the value of d, resulted from the selected pivoting element, is zero or close to zero, then the given matrix AO must be singular.
Anyway the derived algorithm for matrix inversion is very pratical compare to some other classical approaches, such as Gauss-Jordan Elimination. I hope our FEX users will like to use it.
I am not majored in mathematics. I am retired now and soon approach to 80. I need some help regarding the option for automatic pivoting.
The correction works for the reported problem, but, there is more. For instance, for p=[3e-04 1e-09 1e-09 5e-08 2e+04] results are wrong and in addition multiplicities appear as complex numbers.
Reply to Bor Plestenjak:
Thank you for your pinpoint the bug.
It should be OK if the file is correct as follows: In line 15, change
from
" if npo == 1, return, end;"
to
" if npo == 0, return, end;"
I will update this file later. But please try again now!
There seems to be a bug when it comes to linear polynomials (which of course can not have multiple zeros). poly_roots([1 2]) does not return a zero and poly_roots([1 1 0 0 0]) returns just the triple zero 0.
Thank you for all of your comments, John and Jan. Afterall I have decided to continue instead of giving up. Here, I want to say thank again to the encouragement from Jan.
The presentation in File Exchange is now updated to include two files: inv0.m and inv0det_pvt.m.
(1) AI = inv0(AO) -- Find AI(inverse) directly from AO(given matrix) without any pivoting. The code is only 4 statement lines total.
(2) [AI,det] = inv0det_pvt(AO,pvt) -- Find both AI and det(determinant) from AO, with or without pivoting. The "write-in" is used for the pivoting (pvt = 2). It is only good for AO of small order and not practical at all. However it may resolve the problem when the divisor 'd' in the k-loop is zero or close to zero. Such as d=AO(1,1)=0.
The option of automatic pivoting with the preset critique is the way to proceed. I think the critique for pivoting is to select the pivoting element such that the divisor d is the absolute maximum among all computed d's. If the value of d, resulted from the selected pivoting element, is zero or close to zero, then the given matrix AO must be singular.
Anyway the derived algorithm for matrix inversion is very pratical compare to some other classical approaches, such as Gauss-Jordan Elimination. I hope our FEX users will like to use it.
I am not majored in mathematics. I am retired now and soon approach to 80. I need some help regarding the option for automatic pivoting.
Comment only