# [Scilab-users] problem with matrix of polynomials

8 messages
Open this post in threaded view
|

## [Scilab-users] problem with matrix of polynomials

 Dear all, I wonder if the inversion of a matrix of polynomials is working properly. I have the following code from the direct resolution of a two-section constant resistance electrical network:: // Data R1 = 600; R2 = 600; R3 = 600; R4 = 600; R5 = 10070; R6 = 2774; R7 = 36; R8 = 130; R9 = 600; L1 = 2.65; L2 = 0.0318; L3 = 0.00669; L4 = 0.00795; C1 = 18.4e-9; C2 = 22.1e-9; C3 = 7.3e-6; C4 = 88.4e-9; // Series and parallel impedances Z1 = R7 + L3*%s + 1/(C3*%s); Z2 = 1/(1/R5 + 1/(L1*%s) + C1*%s); Z3 = R8 + L4*%s + 1/(C4*%s); Z4 = 1/(1/R6 + 1/(L2*%s) + C2*%s); // System matrix (node potentials) A = [ 1/R1+1/Z1+1/R2, -1/R2,                 0,     0; ...      -1/R2,            1/R2+1/Z2+1/R3+1/Z4, -1/R3, -1/Z4; ...       0,              -1/R3,                 1/R3+1/Z3+1/R4, -1/R4; ...       0,              -1/Z4,                -1/R4,            1/R4+1/Z4+1/R9]; The code goes further, but the point is I need to compute inv(A). When multiplied by A, i.e. A*inv(A) we should get eye(A). Of course, since both A and inv(A) are rationals, we don't get exactly eye(A), but we should anyway get rationals very close to exact num/den cancellation on the diagonal, and orders of magnitude less outside the diagonal. However, it may be a bit tricky to recognize these things in a cluttered output such s the console's, so I attempted replacing values  horner(A*inv(A), %i*2*%pi*1000) I get the following   -894.25452 - 52.256219i   22.488572 - 53.503625i   13.243845 - 53.503625i   13.247333 - 53.503625i    55.248428 - 45.072952i  -905.16326 + 55.026768i  -882.83336 - 2.1937054i  -882.84581 - 2.1939413i    22.907079 - 0.7699342i   877.83532 - 1.5398684i   855.90863 + 56.312642i   878.89745 - 1.5491077i   -64.456283 + 100.49087i  -877.84725 + 1.5398684i  -878.89717 + 1.5491052i  -901.87337 + 59.410849i It doesn't resemble eye(A) at all. Now I tried replacing before inverting: B = horner(inv(A), %i*2*%pi*1000); C = inv(B); Now I compute B*C The result is much closer to eye(A):    1.                       2.845D-16                0.            0.                      -7.216D-16 + 1.332D-15i   1. - 3.886D-16i          1.221D-15i    1.159D-15 + 2.220D-16i    1.249D-16                0.                       1.           -1.366D-16i              -2.776D-16                3.775D-15 - 1.665D-16i   0.            1.                    A last experiment: I replace in the inverse: D = horner(inv(A), %i*2*%pi*1000); Then, computing B*D yields this completely insane result:    1.039D+09 + 7.084D+08i   8.026D+09 + 3.899D+09i   8.103D+09 + 3.626D+09i   8.038D+09 + 3.900D+09i    4.565D+09 - 7.139D+09i   1.425D+11 - 2.222D+10i   1.411D+11 - 2.677D+10i   1.426D+11 - 2.231D+10i   -3.751D+08 - 5.719D+08i  -3.656D+08 - 8.904D+09i  -9.535D+08 - 8.891D+09i  -3.712D+08 - 8.916D+09i    5.187D+09 - 6.874D+09i   1.429D+11 - 2.215D+10i   1.415D+11 - 2.671D+10i   1.431D+11 - 2.225D+10i Am I doing anything wrong?   Regards, Federico Miyara _______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
Open this post in threaded view
|

## [Scilab-users] Re : problem with matrix of polynomials

 Hello Federico, cond() and rcond() do not accept rationals, but the determinant of A is very close to 0: --> det(A)  ans  =                                         3.573D-10              ------------------------      0.0000027s² +7.433D-10s³   This likely explains that the inversion is not reliable. Regards Samuel _______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
Open this post in threaded view
|

## Re: Re : problem with matrix of polynomials

 Maybe an inversion made in state-space could me more reliable ? S. Le 09/03/2020 à 13:36, [hidden email] a écrit : > Hello Federico, > > cond() and rcond() do not accept rationals, but the determinant of A is very close to 0: > --> det(A) >   ans  = >                               >            3.573D-10 >     ------------------------ >     0.0000027s² +7.433D-10s³ > > This likely explains that the inversion is not reliable. > > Regards > Samuel > _______________________________________________ > users mailing list > [hidden email] > https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users-- Stéphane Mottelet Ingénieur de recherche EA 4297 Transformations Intégrées de la Matière Renouvelable Département Génie des Procédés Industriels Sorbonne Universités - Université de Technologie de Compiègne CS 60319, 60203 Compiègne cedex Tel : +33(0)344234688 http://www.utc.fr/~mottelet_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
Open this post in threaded view
|

## Re: Re : problem with matrix of polynomials

 initially I thought could it be possible the number is too small (or too large) to be handled  by double? after trying scaling down the problem, I notice that in Scilab 6.0.2 onwards:--> 2e-20*%i ans  =   0. in Scilab 5.5.2-->2e-20*%i ans  =    2.000D-20i Could this be the issue?thanks.rgds,CL---- On Mon, 09 Mar 2020 21:04:12 +0800 Stéphane Mottelet <[hidden email]> wrote ----Maybe an inversion made in state-space could me more reliable ?S.Le 09/03/2020 à 13:36, [hidden email] a écrit :> Hello Federico,>> cond() and rcond() do not accept rationals, but the determinant of A is very close to 0:> --> det(A)> ans => > 3.573D-10> ------------------------> 0.0000027s² +7.433D-10s³>> This likely explains that the inversion is not reliable.>> Regards> Samuel> _______________________________________________> users mailing list> [hidden email]> https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users-- Stéphane MotteletIngénieur de rechercheEA 4297 Transformations Intégrées de la Matière RenouvelableDépartement Génie des Procédés IndustrielsSorbonne Universités - Université de Technologie de CompiègneCS 60319, 60203 Compiègne cedexTel : +33(0)344234688http://www.utc.fr/~mottelet_______________________________________________users mailing list[hidden email]http://lists.scilab.org/mailman/listinfo/users_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
Open this post in threaded view
|

## Re: &nbsp;&nbsp;&nbsp;&nbsp;Re : problem with matrix of polynomials

 sorry, ignore the previous finding, turn out to be it seems like just a display issue. However, I think my initial thought of the precision might be the issue? some number are too small (first case) and giving error during computation?thxrgds,CL---- On Mon, 09 Mar 2020 22:11:43 +0800 Chin Luh Tan <[hidden email]> wrote ----initially I thought could it be possible the number is too small (or too large) to be handled  by double? after trying scaling down the problem, I notice that in Scilab 6.0.2 onwards:--> 2e-20*%i ans  =   0. in Scilab 5.5.2-->2e-20*%i ans  =    2.000D-20i Could this be the issue?thanks.rgds,CL---- On Mon, 09 Mar 2020 21:04:12 +0800 Stéphane Mottelet <[hidden email]> wrote ----_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users Maybe an inversion made in state-space could me more reliable ?S.Le 09/03/2020 à 13:36, [hidden email] a écrit :> Hello Federico,>> cond() and rcond() do not accept rationals, but the determinant of A is very close to 0:> --> det(A)> ans => > 3.573D-10> ------------------------> 0.0000027s² +7.433D-10s³>> This likely explains that the inversion is not reliable.>> Regards> Samuel> _______________________________________________> users mailing list> [hidden email]> https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users-- Stéphane MotteletIngénieur de rechercheEA 4297 Transformations Intégrées de la Matière RenouvelableDépartement Génie des Procédés IndustrielsSorbonne Universités - Université de Technologie de CompiègneCS 60319, 60203 Compiègne cedexTel : +33(0)344234688http://www.utc.fr/~mottelet_______________________________________________users mailing list[hidden email]http://lists.scilab.org/mailman/listinfo/users_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
Open this post in threaded view
|

## Re: ?==?utf-8?q? ?==?utf-8?q? Re : problem with matrix of polynomials

 In reply to this post by Tan Chin Luh Hello,It seems that this value is below %eps==2.2e-16, so I though that once encoded in double the difference between 2e-20 and 0 is too small to be properly taken into account.Am I missing something?AntoineLe Lundi, Mars 09, 2020 15:11 CET, Chin Luh Tan <[hidden email]> a écrit:  initially I thought could it be possible the number is too small (or too large) to be handled  by double?  after trying scaling down the problem, I notice that in Scilab 6.0.2 onwards: --> 2e-20*%ians  =    0.  in Scilab 5.5.2 -->2e-20*%ians  =     2.000D-20i   Could this be the issue? thanks. rgds,CL ---- On Mon, 09 Mar 2020 21:04:12 +0800 Stéphane Mottelet <[hidden email]> wrote ---- Maybe an inversion made in state-space could me more reliable ?S.Le 09/03/2020 à 13:36, [hidden email] a écrit :> Hello Federico,>> cond() and rcond() do not accept rationals, but the determinant of A is very close to 0:> --> det(A)> ans =>> 3.573D-10> ------------------------> 0.0000027s² +7.433D-10s³>> This likely explains that the inversion is not reliable.>> Regards> Samuel> _______________________________________________> users mailing list> [hidden email]> https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users--Stéphane MotteletIngénieur de rechercheEA 4297 Transformations Intégrées de la Matière RenouvelableDépartement Génie des Procédés IndustrielsSorbonne Universités - Université de Technologie de CompiègneCS 60319, 60203 Compiègne cedexTel : +33(0)344234688http://www.utc.fr/~mottelet_______________________________________________users mailing list[hidden email]http://lists.scilab.org/mailman/listinfo/users   _______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
 In reply to this post by Samuel GOUGEON Samuel, Thanks for your answer. Simply being the determinant small doesn't seem a problem when using floating point. Probably the bigger problem is the presence, somewhere, of very dissimile orders of magnitude. One curious detail in my example is that the solution of the system when adding these lines (input vector and multiplication by the inverse matrix):   V1 = 1; B = [V1/R1, V1/Z2, 0, 0]'; phi = inv(A) * B; the denominator of the resulting node potentials in vector phi should be of degree 8 but they happen to be just scalars (degree 0) as if the coefficients were below %eps Regards, Federico On 09/03/2020 09:36, [hidden email] wrote: ```Hello Federico, cond() and rcond() do not accept rationals, but the determinant of A is very close to 0: --> det(A) ans = 3.573D-10 ------------------------ 0.0000027s² +7.433D-10s³ This likely explains that the inversion is not reliable. Regards Samuel _______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users ``` _______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users