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.
emanuel.c...@gmail.com schrieb am Montag, 10. Oktober 2022 um 19:03:39 UTC+2: > 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/41081428-8b67-46b7-9d20-e336afded913n%40googlegroups.com.
