Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

Unfortunately I cannot find a solution to the problem, but I try to give some information, so that some of you who know better how sagemath works internally may "fix" it (maybe it's a bug in sagemath ?).

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