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

asked 2022-10-07 11:25:01 +0200

Albert_Zweistein
11 ●1 ●2 ●4

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 flag offensive close merge delete

Comments

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.

Max Alekseyev ( 2022-10-07 19:02:51 +0200 )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).

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

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2022-10-09 17:23:34 +0200

Emmanuel Charpentier
7057 ●6 ●46 ●137

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,

edit flag offensive delete link

add a comment

answered 2022-10-08 16:49:16 +0200

Albert_Zweistein
11 ●1 ●2 ●4

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:

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:

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

edit flag offensive delete link

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

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

Comments

2 Answers

Your Answer

Question Tools

Stats

Related questions

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

Comments

2 Answers

Your Answer

Question Tools

Stats

Related questions

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