# Solving a system of linear equations with complex numbers yields false solution

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= 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.

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]


Does anyone know why it yields false results in the first case (taking w*I instead of p) ?

edit retag close merge delete

If you deal with linear equations, it's worth to employ linear algebra machinery (matrices/vectors/etc) rather than symbolic solve function. It'll give you more control over what's going on behind the scene.

( 2022-10-07 19:02:51 +0100 )edit

Maybe you are right. But solving a system of linear equations with constants (real or complex) and variables is a standard task. And since sagemath allows me to do it also without matrices as I did, I would expect to get correct results. Maybe there is some documentation about sagemath that I haven't read and that could explain the behaviour (and in fact I have read only parts so far).

( 2022-10-07 19:42:13 +0100 )edit

Sort by ยป oldest newest most voted

The problem lies in Maxima : Consider the case where I is a variable. Solving it in Maxima :

/* Original system, keeping I as a *variable* */
(%i31) Sys:[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]$(%i32) IVars:[I1, I2, I3, I4]$
/* Solution of the original system */
(%i33) Sol: solve(Sys, IVars)$/* Solution check : */ (%i34) map(lambda([x], is(ratsimp(subst(Sol, x)))), Sys); (%o34) [true, true, true, true]  What happents when I is replaced by the imaginary unit constant %i in Maxima : (%i35) Sys0:map(lambda([x], subst(I=%i, x)), Sys)$
/* Solution */
(%i36) Sol0: solve(Sys0, IVars)\$
/* Checking */
(%i37) map(lambda([x], is(ratsimp(subst(Sol0, x)))), Sys0);
(%o37)                    [true, false, false, true]
(%i38)


The transformed system gets a wrong solution. Somehow, Maxima errs in computations involving the imaginary unit %i as a constant.

This is a bug. Serious, IMHO...

Note : Mathematica does not present this problem :

In Mathematica, I denotes the imaginary unit by default :

In[28]:= I^2
Out[28]= -1

In[29]:= Sys = {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};
In[30]:= IVars = {I1, I2, I3, I4}
Out[30]= {I1, I2, I3, I4}
In[31]:= Sol = Solve[Sys, IVars];


Check (the first element of )this solution :

In[32]:= FullSimplify[#/.Sol[[1]]]&/@Sys
Out[32]= {True, True, True, True}


BTW : this solution is unique:

In[33]:= Length[Sol]
Out[33]= 1


HTH,

more

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

more