[Scilab-users] problem with matrix of polynomials

[Federico Miyara]

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20200309/647f60f4/attachment.htm>