Solving a system of linear equations with complex numbers yields false solution

[Florian Königstein]

I have equations for the analysis of an electric circuit. They contain at several places the term w*I (I is the imaginary unit). I get a solution for the currents I1, I2, I3, I4. Then I substitute the solution into the equations, but I see that two of the four equations are not fulfilled:

var('L1 L2 L3 L4 C1 C2 C3 C4 M23 I1 I2 I3 I4 U w p RL')

lsg = solve([I1*(w*I*L1 + 1/(w*I*C1)) + I2*(-1/(w*I*C1) ) == U, \
             I2*(w*I*L2 + 1/(w*I*C1) + 1/(w*I*C2)) + I3*(w*I*M23 - 1/(w*I*C2)) + I1*(- 1/(w*I*C1)) == 0, \
             I3*(1/(w*I*C3) + w*I*L3 + 1/(w*I*C2)) + I2*(w*I*M23 - 1/(w*I*C2)) - I4/(w*I*C3) == 0, \
             I4/(w*I*C3) + I4*RL - I3/(w*I*C3) == 0], \
             [I1, I2, I3, I4])

param = [w==1, U==1, C1==1, C2==1, C3==1, C4==1, L1==1, L2==1, L3==1, L4==1, RL==1, M23==0]
I1 = I1.subs(lsg).subs(param)
I2 = I2.subs(lsg).subs(param)
I3 = I3.subs(lsg).subs(param)
I4 = I4.subs(lsg).subs(param)

eqn = [I1*(w*I*L1 + 1/(w*I*C1)) + I2*(-1/(w*I*C1) ) == U, \
       I2*(w*I*L2 + 1/(w*I*C1) + 1/(w*I*C2)) + I3*(w*I*M23 - 1/(w*I*C2)) + I1*(- 1/(w*I*C1)) == 0, \
       I3*(1/(w*I*C3) + w*I*L3 + 1/(w*I*C2)) + I2*(w*I*M23 - 1/(w*I*C2)) - I4/(w*I*C3) == 0, \
       I4/(w*I*C3) + I4*RL - I3/(w*I*C3) == 0]

print("I1=", I1)
print("I2=", I2)
print("I3=", I3)
print("I4=", I4)
[eq.subs(param) for eq in eqn]

Output (unexpected):

I1= 1
I2= -I
I3= -2
I4= I - 1
[1 == 1, (-I - 1) == 0, I == 0, 0 == 0]

I have copied and pasted the equations from within the solve() command into the eqn = ... statement. So they are guaranteed to be equal.

Because I have set M23=0 in the parameters, I could remove the terms with M23 in the equations. Then I get the correct solution. I don't understand why not when M23 is present.

Another way to get the correct solution is replacing w*I with p in the equations and setting p=I in the parameters:

var('L1 L2 L3 L4 C1 C2 C3 C4 M23 I1 I2 I3 I4 U w p RL')

lsg = solve([I1*(p*L1 + 1/(p*C1)) + I2*(-1/(p*C1) ) == U, \
             I2*(p*L2 + 1/(p*C1) + 1/(p*C2)) + I3*(p*M23 - 1/(p*C2)) + I1*(- 1/(p*C1)) == 0, \
             I3*(1/(p*C3) + p*L3 + 1/(p*C2)) + I2*(p*M23 - 1/(p*C2)) - I4/(p*C3) == 0, \
             I4/(p*C3) + I4*RL - I3/(p*C3) == 0], \
             [I1, I2, I3, I4])

param = [p==I, U==1, C1==1, C2==1, C3==1, C4==1, L1==1, L2==1, L3==1, L4==1, RL==1, M23==0]
I1 = I1.subs(lsg).subs(param)
I2 = I2.subs(lsg).subs(param)
I3 = I3.subs(lsg).subs(param)
I4 = I4.subs(lsg).subs(param)

eqn = [I1*(p*L1 + 1/(p*C1)) + I2*(-1/(p*C1) ) == U, \
       I2*(p*L2 + 1/(p*C1) + 1/(p*C2)) + I3*(p*M23 - 1/(p*C2)) + I1*(- 1/(p*C1)) == 0, \
       I3*(1/(p*C3) + p*L3 + 1/(p*C2)) + I2*(p*M23 - 1/(p*C2)) - I4/(p*C3) == 0, \
       I4/(p*C3) + I4*RL - I3/(p*C3) == 0]

print("I1=", I1)
print("I2=", I2)
print("I3=", I3)
print("I4=", I4)
[eq.subs(param) for eq in eqn]

Output as expected:

I1= 1
I2= -I
I3= -I - 1
I4= -1
[1 == 1, 0 == 0, 0 == 0, 0 == 0]

When not inserting the numeric parameters, I get the solution in dependence of p (or w), U, C1, C2, L1, L2, M12 and RL. The common denominator of the currents I1, I2, I3, I4 should be the determinant of the matrix (up to a constant factor) if you would write the system in matrix-times-vector form. The denominator is correct both if I use p in the system and if I use I*w instead of p. But the numerators are not correct, at least for I4.

I would like the future of CAS to be in open source like sagemath, but since I used Maple in the past, I begin with some code and results from Maple:

Maple code:
    eqn1 := [I1*(p*L1 + 1/(p*C1)) + I2*(-1/(p*C1) ) = U,
            I2*(p*L2 + 1/(p*C1) + 1/(p*C2)) + I3*(p*M23 - 1/(p*C2)) + I1*(- 1/(p*C1)) = 0,
            I3*(1/(p*C3) + p*L3 + 1/(p*C2)) + I2*(p*M23 - 1/(p*C2)) - I4/(p*C3) = 0,
            I4/(p*C3) + I4*RL - I3/(p*C3) = 0]:
    eqn2 := [seq(subs(p=I*w, eqn1[i]), i=1..4)]:
    lsg1 := solve(eqn1, [I1, I2, I3, I4])[1]:
    lsg2 := solve(eqn2, [I1, I2, I3, I4])[1]:
    print(simplify((subs(p=I*w, subs(lsg1, I4)) - subs(lsg2, I4))));   # expected to be zero

Output is zero as expected.

Maple code:
    print(subs(lsg2, I4));

Output:

-(C2*M23*w^2+1)*U/(-(2*I)*M23*w+(2*I)*C1*L1*M23*w^3-I*C2*M23^2*w^3-RL+C1*L1*RL*w^2+C2*L1*RL*w^2+C3*L1*RL*w^2+C2*L2*RL*w^2+C3*L2*RL*w^2+C3*L3*RL*w^2+2*C3*M23*RL*w^2-I*L1*w-I*L2*w-I*L3*w-C1*C2*L1*L2*RL*w^4-C1*C3*L1*L2*RL*w^4-C1*C3*L1*L3*RL*w^4-C2*C3*L1*L3*RL*w^4-C2*C3*L2*L3*RL*w^4-2*C1*C3*L1*M23*RL*w^4+I*C1*C2*L1*M23^2*w^5-C1*C2*C3*L1*M23^2*RL*w^6-I*C1*C2*L1*L2*L3*w^5+C1*C2*C3*L1*L2*L3*RL*w^6+C2*C3*M23^2*RL*w^4+I*C1*L1*L2*w^3+I*C1*L1*L3*w^3+I*C2*L1*L3*w^3+I*C2*L2*L3*w^3)

Maple code:
    print(numer(subs(lsg2, I4)));       # numerator of I4 from lsg2

Output:

-(C2*M23*w^2+1)*U

Maple code:
    print(denom(subs(lsg2, I4)));       # denominator of I4 from lsg2

Output:

-(2*I)*M23*w+(2*I)*C1*L1*M23*w^3-I*C2*M23^2*w^3-RL+C1*L1*RL*w^2+C2*L1*RL*w^2+C3*L1*RL*w^2+C2*L2*RL*w^2+C3*L2*RL*w^2+C3*L3*RL*w^2+2*C3*M23*RL*w^2-I*L1*w-I*L2*w-I*L3*w-C1*C2*L1*L2*RL*w^4-C1*C3*L1*L2*RL*w^4-C1*C3*L1*L3*RL*w^4-C2*C3*L1*L3*RL*w^4-C2*C3*L2*L3*RL*w^4-2*C1*C3*L1*M23*RL*w^4+I*C1*C2*L1*M23^2*w^5-C1*C2*C3*L1*M23^2*RL*w^6-I*C1*C2*L1*L2*L3*w^5+C1*C2*C3*L1*L2*L3*RL*w^6+C2*C3*M23^2*RL*w^4+I*C1*L1*L2*w^3+I*C1*L1*L3*w^3+I*C2*L1*L3*w^3+I*C2*L2*L3*w^3

Sagemath code:
var('L1 L2 L3 L4 C1 C2 C3 C4 M23 I1 I2 I3 I4 U w p RL')

eqn1 = [I1*(p*L1 + 1/(p*C1)) + I2*(-1/(p*C1) ) == U, \
        I2*(p*L2 + 1/(p*C1) + 1/(p*C2)) + I3*(p*M23 - 1/(p*C2)) + I1*(- 1/(p*C1)) == 0, \
        I3*(1/(p*C3) + p*L3 + 1/(p*C2)) + I2*(p*M23 - 1/(p*C2)) - I4/(p*C3) == 0, \
        I4/(p*C3) + I4*RL - I3/(p*C3) == 0]

eqn2 = [eq.subs(p=I*w) for eq in eqn1]

lsg1 = solve(eqn1, [I1, I2, I3, I4])
lsg2 = solve(eqn2, [I1, I2, I3, I4])

#print((I4.subs(lsg1).subs(p=I*w) - I4.subs(lsg2)))   # expected to be zero but yields long nonzero expression
print(I4.subs(lsg1).subs(p=I*w))

Output:

-(C2*M23*U*w^2 + U)/((C2*C3*L2*L3*RL - C2*C3*M23^2*RL)*C1*L1*w^6 - I*(C2*L2*L3 - C2*M23^2)*C1*L1*w^5 - (C2*C3*L2*L3*RL - C2*C3*M23^2*RL + (C2*C3*L3*RL + (C3*L3*RL + 2*C3*M23*RL + (C2*RL + C3*RL)*L2)*C1)*L1)*w^4 + I*(C2*L2*L3 - C2*M23^2 + (C1*(L2 + L3 + 2*M23) + C2*L3)*L1)*w^3 + (C3*L3*RL + 2*C3*M23*RL + (C1*RL + C2*RL + C3*RL)*L1 + (C2*RL + C3*RL)*L2)*w^2 - I*(L1 + L2 + L3 + 2*M23)*w - RL)

Sagemath code:
    print(I4.subs(lsg2))

Output:

(I*C2*C3*M23*RL*U*w^3 + C2*M23*U*w^2 + I*C3*RL*U*w + U)/(-I*(C2*C3*L2*L3*RL - C2*C3*M23^2*RL)*C1*L1*w^6 - (C2*L2*L3 - C2*M23^2)*C1*L1*w^5 + (I*C2*C3*L2*L3*RL - I*C2*C3*M23^2*RL + (I*C2*C3*L3*RL + I*(C3*L3*RL + 2*C3*M23*RL + (C2*RL + C3*RL)*L2)*C1)*L1)*w^4 + (C2*L2*L3 - C2*M23^2 + (C1*(L2 + L3 + 2*M23) + C2*L3)*L1)*w^3 + (-I*C3*L3*RL - 2*I*C3*M23*RL + (-I*C1*RL - I*C2*RL - I*C3*RL)*L1 - I*(C2*RL + C3*RL)*L2)*w^2 - (L1 + L2 + L3 + 2*M23)*w + I*RL)

Sagemath code:
print((denominator(I4.subs(lsg1)).subs(p=I*w) / denominator(I4.subs(lsg2))).simplify_full())

Output:

I

Sagemath code:
print((denominator(I4.subs(lsg1)).subs(p=I*w) - I*denominator(I4.subs(lsg2))).simplify_full())  # yields zero since previous quotiont is I

Output:

0

Sagemath code:
print((numerator(I4.subs(lsg1)).subs(p=I*w) - I*numerator(I4.subs(lsg2))).simplify_full())  # expected to be zero

Output:

-C2*C3*M23*RL*U*w^3 + (I + 1)*C2*M23*U*w^2 - C3*RL*U*w + (I + 1)*U

Sagemath code:
print(numerator(I4.subs(lsg1)).subs(p=I*w))

Output:

C2*M23*U*w^2 + U

Sagemath code:
print(numerator(I4.subs(lsg2)))

Output:

-I*C2*C3*M23*RL*U*w^3 - C2*M23*U*w^2 - I*C3*RL*U*w - U

For me it looks that it could be a bug in sagemath. If not, I would like to know why sagemath behaves like this.

[David Joyner]

I think you essentially answered your own question. IMHO, the symbolic
solver implicitly assumes a symbolic parameter is non-zero when
solving for the other variables. It looks like your equations are
linear in the variables to be solved for and so the solver is using
row reduction at some point., which of course involves divisions by
(presumably non-zero) expressions in the parameters.

> --
> You received this message because you are subscribed to the Google Groups "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/8d56151f-5295-455c-a86e-92ae1a04cf88n%40googlegroups.com.

[Florian Königstein]

David Joyner schrieb am Sonntag, 9. Oktober 2022 um 10:23:51 UTC+2:

IMHO, the symbolic solver implicitly assumes a symbolic parameter is non-zero when
solving for the other variables.

It looks like your equations are
linear in the variables to be solved for and so the solver is using
row reduction at some point., which of course involves divisions by
(presumably non-zero) expressions in the parameters.
 

Yes, the equations are linear in I1, I2, I3, I4. In matrix-vector form they are:

var('L1 L2 L3 C1 C2 C3 M23 I1 I2 I3 I4 U w p RL')
A1 = Matrix(
[
    [p*L1 + 1/(p*C1) , -1/(p*C1)                                     , 0                                        , 0            ],
    [-1/(p*C1)            , p*L2 + 1/(p*C1) + 1/(p*C2)      ,  p*M23 - 1/(p*C2)          , 0            ],
    [0                          , p*M23 - 1/(p*C2)                       , 1/(p*C3) + p*L3 + 1/(p*C2) , -1/(p*C3)    ],
    [0                          , 0                                                    , -1/(p*C3)                          , 1/(p*C3) + RL]
])
RightSide = vector([U, 0, 0, 0])
A2 = A1.subs(p=w*I)

Currents1 = A1.solve_right(RightSide)
Currents2 = A2.solve_right(RightSide)

# verify:

print((A1*Currents1).simplify_full())

Output as expected : (U, 0, 0, 0)

print((A2*Currents2).subs(param))

Output as expected : (U, 0, 0, 0)

So, with the matrix-vector formalism is works correctly both when I use p and when I use w*I instead.

But IMHO when using the solve-command from my first post with using w*I instead of p it should also work

correctly - when I do not insert numerical parameters for L1, L2, L3, C1, C2, C3, M23, p and w and when I perform

a .simplify_all() on the result.

Assume I use the solve()-command with w*I instead of p and the solver uses the Gauss-Seidel-algorithm to solve for I1, I2, I3, I4.

This corresponds to the matrix A2 in this post. Assume the solver chooses first the "pivot"-element A2[2,1] . This element is I*M23*w + I/(C2*w) .

The solver doesn't know numeric values for M23, w and C2 - but it assumes that A2[2,1] is nonzero. The solver continues until having a

solution. If I insert then numeric values, then it could in fact be that A2[2,1] gets zero and so you could explain wrong results. But I

can insert other numeric values with A2[2,1] nonzero and then the pivot element was choosen correctly.

When I don't insert numeric values and perform a simplify_full() on I1, I2, I3, I4, I get rational functions for them: The common denominator

is equal to the determinant of A2 (maybe up to a constant factor that is independent from L1, L2, L3, C1, C2, C3, M23, RL, w (the variables inside of A2)).

The numerators of I1, I2, I3, I4 can e.g. be found by Cramer's rule.

IMHO using the solve-command with w*I instead of p and applying a simplify_all() should yield the same results as if you determined the

rational functions I1, I2, I3, I4 with Cramer's rule. The only case when it can and must fail is when the determinant of A2 is zero after inserting

numeric values for L1, L2, L3, C1, C2, C3, M23, RL and w .

A nonzero determinant is e.g. generated by:

param = [w==1, p==I, C1==2, C2==3, C3==4, C4==5, L1==6, L2==7, L3==8, L4==9, RL==10, M23==1]

I come back to the case that A2[2,1]==0 for special numeric values, but I assume that the determinant of A2 is nonzero for these values.

If the solver has choosen A2[2,1] as first "pivot" element, the division by zero ( 1 / A[2,1] ) should be cancelled out by a multiplication

with A[2,1] by the simplify_full() command.

[Dima Pasechnik]

On Sun, Oct 9, 2022 at 9:23 AM David Joyner <wdjo...@gmail.com> wrote:
>
> On Sun, Oct 9, 2022 at 3:25 AM Florian Königstein <niets...@gmail.com> wrote:
> >
> > I have equations for the analysis of an electric circuit. They contain at several places the term w*I (I is the imaginary unit). I get a solution for the currents I1, I2, I3, I4. Then I substitute the solution into the equations, but I see that two of the four equations are not fulfilled:
> >
> > var('L1 L2 L3 L4 C1 C2 C3 C4 M23 I1 I2 I3 I4 U w p RL')
> >
> > lsg = solve([I1*(w*I*L1 + 1/(w*I*C1)) + I2*(-1/(w*I*C1) ) == U, \
> > I2*(w*I*L2 + 1/(w*I*C1) + 1/(w*I*C2)) + I3*(w*I*M23 - 1/(w*I*C2)) + I1*(- 1/(w*I*C1)) == 0, \
> > I3*(1/(w*I*C3) + w*I*L3 + 1/(w*I*C2)) + I2*(w*I*M23 - 1/(w*I*C2)) - I4/(w*I*C3) == 0, \
> > I4/(w*I*C3) + I4*RL - I3/(w*I*C3) == 0], \
> > [I1, I2, I3, I4])

This is a homogeneous 4x4 linear system, it will only have nonzero
solutions where the determinant of its matrix
vanishes. If you'd like to analyse the determinant, best is to work
with polynomials, not symbolic variables.

ii=QQ['ii'].0
C.<ii>=NumberField(ii^2+1)
K.<w,L1,C1,L2,C2,M23,C3,L3,RL>=C[] # replacing I by ii, I is too
special in Sage...
# I1 I2 I3 I4
A=matrix(K.fraction_field(), [
[w*ii*L1+1/(w*ii*C1),-1/(w*ii*C1), 0,
0 ],
[-1/(w*ii*C1), w*ii*L2+1/(w*ii*C1)+1/(w*ii*C2), w*ii*M23 -
1/(w*ii*C2), 0 ],
[0, w*ii*M23-1/(w*ii*C2),
1/(w*ii*C3)+w*ii*L3+1/(w*ii*C2), -1/(w*ii*C3) ],
[0, 0, -1/(w*ii*C3),
1/(w*ii*C3) + RL]
])


sage: load('x.sage')
sage: A.det().factor()
(ii) * w^-3 * C3^-1 * C2^-1 * C1^-1 * (w^6*L1*C1*C2*M23^2*C3*RL -
w^6*L1*C1*L2*C2*C3*L3*RL + (-ii)*w^5*L1*C1*C2*M23^2 +
ii*w^5*L1*C1*L2*C2*L3 + w^4*L1*C1*L2*C2*RL + w^4*L1*C1*L2*C3*RL +
2*w^4*L1*C1*M23*C3*RL - w^4*C2*M23^2*C3*RL + w^4*L1*C1*C3*L3*RL +
w^4*L1*C2*C3*L3*RL + w^4*L2*C2*C3*L3*RL + (-ii)*w^3*L1*C1*L2 +
(-2*ii)*w^3*L1*C1*M23 + ii*w^3*C2*M23^2 + (-ii)*w^3*L1*C1*L3 +
(-ii)*w^3*L1*C2*L3 + (-ii)*w^3*L2*C2*L3 - w^2*L1*C1*RL - w^2*L1*C2*RL
- w^2*L2*C2*RL - w^2*L1*C3*RL - w^2*L2*C3*RL - 2*w^2*M23*C3*RL -
w^2*C3*L3*RL + ii*w*L1 + ii*w*L2 + (2*ii)*w*M23 + ii*w*L3 + RL)

You haven't told us what ranges your parameters have (some constraints
can be seen from the matrix, e.g. C1, C2, C3, and w cannot be 0).
E.g. if they are all reals then det(A)=0 is actually 2 equations, one
for real and one for the imaginary part.

It's also remarkable that A is tridiagonal, and so one can rather
easily write, assuming I1 nonzero (if I1=0 then I2=0, and this is a
much easier problem),
an expression for I2 in terms of I1 (and entries of 1st row of A),
I3 in terms of I1, I2, and the 1st and 2nd row of A (assuming
w*ii*M23 - 1/(w*ii*C2) is not zero - if it is zero, then you have two
indendent blocks in A, again a much
easier problem).
Finally, I4 can be expressed in terms of I1,I2,I3 and entries of A in
two different ways - one using the 3rd row of A, another using 4th row
of A.
These two expressions must obviosuly be equal (which will only happen
if det(A)=0).
So it seems one can gather a lot of info about all the possible
solutions, or even write down a full answer.

Hope this helps,
Dima

> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/CAEQuuAWVd%2Bg7xy74NVA40oYuwHPGL3oys6iz42htAgSMCO%3Db5Q%40mail.gmail.com.

[Emmanuel Charpentier]

I posted about the same problem ; the original post in ask.sagemath.org is clearer. The original problem seems to emane from Maxima.

​

[Dima Pasechnik]

Sorry, I mistook U in the RHS of the 1st equation for 0, and so what I wrote only applies to the case U=0.

For U nonzero, one can write a generic solution, with A being full rank, using the Cramer rule.

Although there could be more solutions with det(A)=0, and U being nonzero, as well.

[Florian Königstein]

Thank you for your efforts.

Yes, it was me who posted the original post https://ask.sagemath.org/question/64344/solving-a-system-of-linear-equations-with-complex-numbers-yields-false-solution/  .

I think  my member name "Albert Zweistein" (allusion to Albert Einstein) is a bit snobbish. Is there a way to change the name in ask.sagemath.org ?

I received your posts while I was writing my post.

Since it seems to be clear now that it's a bug, I only write (as additional information) from where I got the equations.

But for fixing the bug the following information in absolutely not necessary:

For those who have some knowledge of analysis of electrical circuits with complex numbers: In the attachment I have an image
of the circuit from that I have generated the equations with the mesh current analysis. The physical time-dependent values for
voltages and currents are of course real, but you calculate with complex numbers: If the driving voltage is described with the complex
number U, the physical voltage is by definition the real part of U*exp(I*w*t) , where t is the time. Correspondingly e.g. the physical current i1 is
equal to the real part of I1*exp(I*w*t).

Physically the inductances L1, L2, L3 and the mutual inductance M23 are real and positive, but for design purposes you could also
allow them to get negative (then you would replace a negative inductance by a capacitor if you know the frequency of the driving voltage).
The capacitances C1, C2, C3 are also real and normally positive, but for design purposes you could also allow negative values.
The resistance RL is normally positive (or may be zero).
The angular frequency w is real and nonnegative, and so p = I*w is purely imaginary.

[Emmanuel Charpentier]

This problem has been filed in Maxima’s bug report system.

Should a Sage ticket be filed ?

​

[Florian Königstein]

I think it would make sense to file a Sage ticket. So when it will be fixed in Maxima, it can be verified that it (hopefully) will also work in Sage with the updated version of Maxima.

[Emmanuel Charpentier]

Le lundi 10 octobre 2022 à 20:37:57 UTC+2, Florian Königstein a écrit :

I think it would make sense to file a Sage ticket. So when it will be fixed in Maxima, it can be verified that it (hopefully) will also work in Sage with the updated version of Maxima.

Done.

[Florian Königstein]

OK, thank you for your efforts.

[Florian Königstein]

I'd like to point out that I found a much simpler system of equations for that it fails. I have also reported in the thread for maxima: https://sourceforge.net/p/maxima/bugs/4032/

var('I1 I2 I4 I5 C1 w')
assume(I5, "real", I5>0)
assume(C1, "real", C1>0)
assume(w, "real", w>0)

sys = [I*I1*w*C1 == 1, I2*w == I5*w, I*I4*w + I4 == 0]
lsg = solve(sys, [I1,I2,I4])
print(lsg)
[bool(eq.subs(lsg).simplify_full()) for eq in sys]

The output is:

[ [I1 == (-I*w^2 - 2*w + I)/(C1*w^2 - I*C1*w), I2 == I5, I4 == 0] ]

[False, True, True]

The system matrix is already diagonal, so that the system can be solved without any Gauss elimination steps. But the solution for I1 is wrong.

[Emmanuel Charpentier]

FWIW, solve(..., algorithm="sympy") solves this simplified system.

HTH,

​

[dmo...@deductivepress.ca]

Thanks for the simpler example.  I copied it to the trac ticket (#34650).  Please post any further comments to the trac ticket, instead of to this discussion.