This problem has been filed <https://sourceforge.net/p/maxima/bugs/4032/> in Maxima’s bug report system.
Should a Sage ticket be filed ? Le dimanche 9 octobre 2022 à 19:55:07 UTC+2, Florian Königstein a écrit : > 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. > [image: Circuit.png] > dim...@gmail.com schrieb am Sonntag, 9. Oktober 2022 um 18:55:26 UTC+2: > >> 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. >> >> >> On Sunday, October 9, 2022 at 1:15:57 PM UTC+1 Dima Pasechnik wrote: >> >>> 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/e0b32aee-3338-435d-aebc-6c213e33d261n%40googlegroups.com.
