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.
Hi, I have a polynomial with 34 orders, but this function gives me roots of 7x2 matrix. Why is that? I copied my polynomial coefficients as below. Could anyone help me out? Thanks!
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.
Hi, I have a polynomial with 34 orders, but this function gives me roots of 7x2 matrix. Why is that? I copied my polynomial coefficients as below. Could anyone help me out? Thanks!
1 -26.7479290179817 346.353790281681 -2891.79278087362 17494.1298230978 -81695.4463236101 306371.295805817 -947699.866103653 2464717.23435414 -5466000.43064478 10447379.3286537 -17350483.7757060 25192808.9914633 -32131544.3997781 36120714.4912932 -35872219.7954520 31515357.4741505 -24504114.0520695 16854492.5963053 -10241648.1386483 5485377.96325360 -2580928.16918672 1061956.10387350 -379845.324952006 117197.882178068 -30882.3382140184 6860.29267091761 -1262.93321056908 188.263121608582 -21.9982275506094 1.92027369835658 -0.115917764563237 0.00419847428401219 -6.48118467521201e-05
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.
Comment only