I posted <https://groups.google.com/g/sage-support/c/V_vdS2ipgek> about the
same problem ; the original post
<https://ask.sagemath.org/question/64344/solving-a-system-of-linear-equations-with-complex-numbers-yields-false-solution/>
in ask.sagemath.org is clearer. The original problem seems to emane from
Maxima.
Le dimanche 9 octobre 2022 à 14:15:57 UTC+2, dim...@gmail.com a écrit :
> 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
>
>
>
>
>
> > >
> > > 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.
> > >
> >
> > 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.
> >
> > > 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.
> > >
> > > --
> > > 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
> .
> >
> > --
> > 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/CAEQuuAWVd%2Bg7xy74NVA40oYuwHPGL3oys6iz42htAgSMCO%3Db5Q%40mail.gmail.com
> .
>
--
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sage-devel/8396b679-7c85-4fa4-b154-c6004728f353n%40googlegroups.com.
