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= 1I2= -II3= -2I4= 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= 1I2= -II3= -I - 1I4= -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 expressionprint(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+unsubscr...@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.
