[Scilab-users] problem with matrix of polynomials

classic Classic list List threaded Threaded
8 messages Options
fmiyara fmiyara
Reply | Threaded
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
Samuel GOUGEON Samuel GOUGEON
Reply | Threaded
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
mottelet mottelet
Reply | Threaded
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
Tan Chin Luh Tan Chin Luh
Reply | Threaded
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 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



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Tan Chin Luh Tan Chin Luh
Reply | Threaded
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?

thx

rgds,
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 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



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Antoine Monmayrant-2 Antoine Monmayrant-2
Reply | Threaded
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?

Antoine


Le 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*%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 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
 

 
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Tan Chin Luh Tan Chin Luh
Reply | Threaded
Open this post in threaded view
|

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

Hi, 

sorry for this misleading finding, infact I've second message to turn it down.


******************
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?
******************

Back to the possible reason of causing this issue when u substitute a value into the %s  into A*inv(A)

Assume a = A*inv(A), a(1,1) which suppose to be ~ 1 is : 

--> a(1,1)
ans  =

                                                                   
   1.630D-38 + 1.586D-40s + 4.110D-43s^2 + 1.120D-46s^3 + 1.998D-50s^4 + 3.586D-55s^5 
   --------------------------------------------------------------------------- 
                                                                          
                      1.630D-38 + 4.283D-42s + 7.960D-46s^2

the coefficients in numerator and denominator are very small value and with complex number involve, there might be come kind of limitation in the numerical software. 

this could be shown if you subs with smaller number: 

--> horner(A*inv(A), %i*2*%pi*1)
ans  =

   0.9991047 + 0.0595012i   0.0000225 - 5.350D-08i   0.0000132 - 5.350D-08i   0.0000132 - 5.350D-08i
   0.0000552 - 0.1037974i   0.9990938 + 0.0595013i  -0.0008828 - 2.194D-09i  -0.0008828 - 2.194D-09i
   0.0000229 - 7.699D-10i   0.0008778 - 1.540D-09i   1.0008549 + 0.0595013i   0.0008789 - 1.549D-09i
  -0.0000645 + 0.1025666i  -0.0008778 + 1.540D-09i  -0.0008789 + 1.549D-09i   0.9990971 + 0.0595013i

you get the correct answer

This is so far I could go, further details such as how the a double precision real or complex number handled in numerical software are beyond my knowledge.
thanks.

Regards,
Chin Luh




---- On Mon, 09 Mar 2020 22:20:38 +0800 Antoine Monmayrant <[hidden email]> wrote ----

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?

Antoine


Le 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*%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 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
 

  _______________________________________________
users mailing list



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
fmiyara fmiyara
Reply | Threaded
Open this post in threaded view
|

Re: Re : problem with matrix of polynomials

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