sympy solve() gives implicit/incorrect answer

Question

I'm trying to solve an equation system with 16 equations and 16 unknowns using sympy but it doesn't seem to solve it well.

I want to solve the system [K][d]=[f] where [K] is the coefficients matrix, [d] the unknowns and [f] are constants. I know some unknowns "d" and some constants "f", so I have same number for both equations and unknowns, but when I substitute these values into the equations and try to solve it the results for all "dx" include "dx8". I checked the matrix determinant and is positive so I should get a unique answer.

Here is the code:

import sympy as sp
import numpy as np

K = np.array([[560000000.0, 0.0, -480000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0,-80000000.0, 120000000.0, 0.0, -200000000.0, 0.0, 0.0, 0.0, 0.0],
   [0.0, 393333333.3, 120000000.0, -180000000.0, 0.0, 0.0, 0.0, 0.0,80000000.0, -213333333.3, -200000000.0, 0.0, 0.0, 0.0, 0.0, 0.0],
   [-480000000.0, 120000000.0, 1120000000.0, -200000000.0,-480000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0, -160000000.0,200000000.0, 0.0, -200000000.0, 0.0, 0.0],
   [80000000.0, -180000000.0, -200000000.0, 786666666.7, 120000000.0,-180000000.0, 0.0, 0.0, 0.0, 0.0, 200000000.0, -426666666.7,-200000000.0, 0.0, 0.0, 0.0],
   [0.0, 0.0, -480000000.0, 120000000.0, 1120000000.0, -200000000.0,-480000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0, -160000000.0,200000000.0, 0.0, -200000000.0],
   [0.0, 0.0, 80000000.0, -180000000.0, -200000000.0, 786666666.7,120000000.0, -180000000.0, 0.0, 0.0, 0.0, 0.0, 200000000.0,-426666666.7, -200000000.0, 0.0],
   [0.0, 0.0, 0.0, 0.0, -480000000.0, 120000000.0, 560000000.0,-200000000.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -80000000.0, 80000000.0],
   [0.0, 0.0, 0.0, 0.0, 80000000.0, -180000000.0, -200000000.0,393333333.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 120000000.0,-213333333.3],
   [-80000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 560000000.0,-200000000.0, -480000000.0, 120000000.0, 0.0, 0.0, 0.0, 0.0],
   [120000000.0, -213333333.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-200000000.0, 393333333.3, 80000000.0, -180000000.0, 0.0, 0.0, 0.0,0.0],
   [0.0, -200000000.0, -160000000.0, 200000000.0, 0.0, 0.0, 0.0, 0.0,-480000000.0, 80000000.0, 1120000000.0, -200000000.0, -480000000.0,120000000.0, 0.0, 0.0],
   [-200000000.0, 0.0, 200000000.0, -426666666.7, 0.0, 0.0, 0.0, 0.0,120000000.0, -180000000.0, -200000000.0, 786666666.7, 80000000.0,-180000000.0, 0.0, 0.0],
   [0.0, 0.0, 0.0, -200000000.0, -160000000.0, 200000000.0, 0.0, 0.0,0.0, 0.0, -480000000.0, 80000000.0, 1120000000.0, -200000000.0,-480000000.0, 120000000.0],
   [0.0, 0.0, -200000000.0, 0.0, 200000000.0, -426666666.7, 0.0, 0.0,0.0, 0.0, 120000000.0, -180000000.0, -200000000.0, 786666666.7,80000000.0, -180000000.0],
   [0.0, 0.0, 0.0, 0.0, 0.0, -200000000.0, -80000000.0, 120000000.0,0.0, 0.0, 0.0, 0.0, -480000000.0, 80000000.0, 560000000.0, 0.0],
   [0.0, 0.0, 0.0, 0.0, -200000000.0, 0.0, 80000000.0, -213333333.3,0.0, 0.0, 0.0, 0.0, 120000000.0, -180000000.0, 0.0, 393333333.3]])



x = [sp.var('dx'+ str(i+1)) for i in range(8)]
y = [sp.var('dy'+ str(i+1)) for i in range(8)]
fx = [sp.var('fx'+ str(i+1)) for i in range(8)]
fy = [sp.var('fy'+ str(i+1)) for i in range(8)]

xy = list(sum(zip(x, y), ()))
fxy = list(sum(zip(fx, fy), ()))

M = sp.Matrix(K)*sp.Matrix(xy)
Ec = [sp.Eq(M[i], fxy[i]) for i in range(16)]

#known values
d_kwn = [(dy1, 0), (dy2, 0), (dy3, 0), (dy4, 0)]

f_kwn = [(fx5, 0), (fy5, 0), (fx6, 0), (fy6, -3000), (fx7, 0), (fy7, -3000),(fx8, 0), (fy8, 0), (fx1, 0), (fx2, 0), (fx3, 0), (fx4, 0)]

for var in d_kwn:
    for i, eq in enumerate(Ec):
        Ec[i] = eq.subs(var[0], var[1])

for var in f_kwn:
    for i, eq in enumerate(Ec):
        Ec[i] = eq.subs(var[0], var[1])

Sols = sp.solvers.solve(Ec)
sp.Matrix(sorted(Sols.items(), key=str))

And this is the output I'm getting:

{dx1:  dx8−3.54468009860439⋅10−6,
dx2:  dx8−1.8414987360977⋅10−6,
dx3:  dx8−2.11496606381994⋅10−7,
dx4:  dx8+2.05943267588118⋅10−7,
dx5:  dx8−1.24937663359153⋅10−6,
dx6:  dx8−1.55655946713284⋅10−6,
dx7:  dx8−1.08797652070783⋅10−6,
dy5:  −2.10639657360695⋅10−6,
dy6:  −6.26959460018537⋅10−6,
dy7:  −6.32191585665888⋅10−6,
dy8:  −2.7105825114088⋅10−6,
fy1:  439.746516706791,
fy2:  2640.65618690176,
fy3:  2399.44807607611,
fy4:  520.14922031534}

I don't know why I'm not getting a result for dx8. I tried adding more equations because theoretically: dx1 = dx4, dx2 = dx3, dx5 = dx8, dx6 = dx7 and so on. But it gives me and empty list. Any help will be appreciated.

Answer 1

It is convenient to solve system of linear equations in Numpy itself. The type of system you are solving appears in Finite Element Analysis often with boundary conditions. Is it fine if we only use Numpy? If yes, the following code will do the job. We already know which elements of f and d are known we can use Numpy array indexing to solve the reduced set of equations as follows:

import numpy as np

# The NxN Coefficients matrix
K = np.array([[560000000.0, 0.0, -480000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0,-80000000.0, 120000000.0, 0.0, -200000000.0, 0.0, 0.0, 0.0, 0.0],
   [0.0, 393333333.3, 120000000.0, -180000000.0, 0.0, 0.0, 0.0, 0.0,80000000.0, -213333333.3, -200000000.0, 0.0, 0.0, 0.0, 0.0, 0.0],
   [-480000000.0, 120000000.0, 1120000000.0, -200000000.0,-480000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0, -160000000.0,200000000.0, 0.0, -200000000.0, 0.0, 0.0],
   [80000000.0, -180000000.0, -200000000.0, 786666666.7, 120000000.0,-180000000.0, 0.0, 0.0, 0.0, 0.0, 200000000.0, -426666666.7,-200000000.0, 0.0, 0.0, 0.0],
   [0.0, 0.0, -480000000.0, 120000000.0, 1120000000.0, -200000000.0,-480000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0, -160000000.0,200000000.0, 0.0, -200000000.0],
   [0.0, 0.0, 80000000.0, -180000000.0, -200000000.0, 786666666.7,120000000.0, -180000000.0, 0.0, 0.0, 0.0, 0.0, 200000000.0,-426666666.7, -200000000.0, 0.0],
   [0.0, 0.0, 0.0, 0.0, -480000000.0, 120000000.0, 560000000.0,-200000000.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -80000000.0, 80000000.0],
   [0.0, 0.0, 0.0, 0.0, 80000000.0, -180000000.0, -200000000.0,393333333.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 120000000.0,-213333333.3],
   [-80000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 560000000.0,-200000000.0, -480000000.0, 120000000.0, 0.0, 0.0, 0.0, 0.0],
   [120000000.0, -213333333.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-200000000.0, 393333333.3, 80000000.0, -180000000.0, 0.0, 0.0, 0.0,0.0],
   [0.0, -200000000.0, -160000000.0, 200000000.0, 0.0, 0.0, 0.0, 0.0,-480000000.0, 80000000.0, 1120000000.0, -200000000.0, -480000000.0,120000000.0, 0.0, 0.0],
   [-200000000.0, 0.0, 200000000.0, -426666666.7, 0.0, 0.0, 0.0, 0.0,120000000.0, -180000000.0, -200000000.0, 786666666.7, 80000000.0,-180000000.0, 0.0, 0.0],
   [0.0, 0.0, 0.0, -200000000.0, -160000000.0, 200000000.0, 0.0, 0.0,0.0, 0.0, -480000000.0, 80000000.0, 1120000000.0, -200000000.0,-480000000.0, 120000000.0],
   [0.0, 0.0, -200000000.0, 0.0, 200000000.0, -426666666.7, 0.0, 0.0,0.0, 0.0, 120000000.0, -180000000.0, -200000000.0, 786666666.7,80000000.0, -180000000.0],
   [0.0, 0.0, 0.0, 0.0, 0.0, -200000000.0, -80000000.0, 120000000.0,0.0, 0.0, 0.0, 0.0, -480000000.0, 80000000.0, 560000000.0, 0.0],
   [0.0, 0.0, 0.0, 0.0, -200000000.0, 0.0, 80000000.0, -213333333.3,0.0, 0.0, 0.0, 0.0, 120000000.0, -180000000.0, 0.0, 393333333.3]])

# A logical array for indexing
N = K.shape[0] # The number of columns in K
N_2 = int(N/2);

# Prepare the 'f'
fx = np.zeros( N_2 );
fy = np.zeros( N_2 );

fx[ [0,1,2,3,4,5,6,7] ] = np.array([0]*N_2) # Known values of fx
fy[ [4,5,6,7] ] = np.array([0,-3000,-3000,0])

f = np.concatenate( (fx,fy) )

# Solve for the unknown equations only
d = np.zeros( N )
rows = np.array([0,1,2,3,4,5,6,7,12,13,14,15])
rows = rows[:, np.newaxis]
columns = np.array([0,1,2,3,4,5,6,7,12,13,14,15])

d[ columns ] = np.linalg.solve( K[ rows, columns ], f[ columns ] )

# Calculate unknown f values
f[ [8,9,10,11] ] = K[ [8,9,10,11], [8,9,10,11] ]*d[[8,9,10,11]]

Answer 2

If you need to use Sympy, then the following may work. First we can solve the reduced system of equations only for unknown d values. Then once we know all d values we can calculate the unknown f values by doing [K][d]=[f] for only the unknown f equation numbers (not implemented in the code below).

import sympy as sp
import numpy as np

K = np.array([[560000000.0, 0.0, -480000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0,-80000000.0, 120000000.0, 0.0, -200000000.0, 0.0, 0.0, 0.0, 0.0],
   [0.0, 393333333.3, 120000000.0, -180000000.0, 0.0, 0.0, 0.0, 0.0,80000000.0, -213333333.3, -200000000.0, 0.0, 0.0, 0.0, 0.0, 0.0],
   [-480000000.0, 120000000.0, 1120000000.0, -200000000.0,-480000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0, -160000000.0,200000000.0, 0.0, -200000000.0, 0.0, 0.0],
   [80000000.0, -180000000.0, -200000000.0, 786666666.7, 120000000.0,-180000000.0, 0.0, 0.0, 0.0, 0.0, 200000000.0, -426666666.7,-200000000.0, 0.0, 0.0, 0.0],
   [0.0, 0.0, -480000000.0, 120000000.0, 1120000000.0, -200000000.0,-480000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0, -160000000.0,200000000.0, 0.0, -200000000.0],
   [0.0, 0.0, 80000000.0, -180000000.0, -200000000.0, 786666666.7,120000000.0, -180000000.0, 0.0, 0.0, 0.0, 0.0, 200000000.0,-426666666.7, -200000000.0, 0.0],
   [0.0, 0.0, 0.0, 0.0, -480000000.0, 120000000.0, 560000000.0,-200000000.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -80000000.0, 80000000.0],
   [0.0, 0.0, 0.0, 0.0, 80000000.0, -180000000.0, -200000000.0,393333333.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 120000000.0,-213333333.3],
   [-80000000.0, 80000000.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 560000000.0,-200000000.0, -480000000.0, 120000000.0, 0.0, 0.0, 0.0, 0.0],
   [120000000.0, -213333333.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-200000000.0, 393333333.3, 80000000.0, -180000000.0, 0.0, 0.0, 0.0,0.0],
   [0.0, -200000000.0, -160000000.0, 200000000.0, 0.0, 0.0, 0.0, 0.0,-480000000.0, 80000000.0, 1120000000.0, -200000000.0, -480000000.0,120000000.0, 0.0, 0.0],
   [-200000000.0, 0.0, 200000000.0, -426666666.7, 0.0, 0.0, 0.0, 0.0,120000000.0, -180000000.0, -200000000.0, 786666666.7, 80000000.0,-180000000.0, 0.0, 0.0],
   [0.0, 0.0, 0.0, -200000000.0, -160000000.0, 200000000.0, 0.0, 0.0,0.0, 0.0, -480000000.0, 80000000.0, 1120000000.0, -200000000.0,-480000000.0, 120000000.0],
   [0.0, 0.0, -200000000.0, 0.0, 200000000.0, -426666666.7, 0.0, 0.0,0.0, 0.0, 120000000.0, -180000000.0, -200000000.0, 786666666.7,80000000.0, -180000000.0],
   [0.0, 0.0, 0.0, 0.0, 0.0, -200000000.0, -80000000.0, 120000000.0,0.0, 0.0, 0.0, 0.0, -480000000.0, 80000000.0, 560000000.0, 0.0],
   [0.0, 0.0, 0.0, 0.0, -200000000.0, 0.0, 80000000.0, -213333333.3,0.0, 0.0, 0.0, 0.0, 120000000.0, -180000000.0, 0.0, 393333333.3]])



x = [sp.var('dx'+ str(i+1)) for i in range(8)]
y = [sp.var('dy'+ str(i+1)) for i in range(8)]
fx = [sp.var('fx'+ str(i+1)) for i in range(8)]
fy = [sp.var('fy'+ str(i+1)) for i in range(8)]

xy = list(sum(zip(x, y), ()))
fxy = list(sum(zip(fx, fy), ()))

M = sp.Matrix(K)*sp.Matrix(xy)
Ec = [sp.Eq(M[i], fxy[i]) for i in range(16)]

#known values
d_kwn = [(dy1, 0), (dy2, 0), (dy3, 0), (dy4, 0)]

f_kwn = [(fx5, 0), (fy5, 0), (fx6, 0), (fy6, -3000), (fx7, 0), (fy7, -3000),(fx8, 0), (fy8, 0), (fx1, 0), (fx2, 0), (fx3, 0), (fx4, 0)]

for var in d_kwn:
    for i, eq in enumerate(Ec):
        Ec[i] = eq.subs(var[0], var[1])

for var in f_kwn:
    for i, eq in enumerate(Ec):
        Ec[i] = eq.subs(var[0], var[1])

Ec_part = []
for i in [0,2,4,6,8,9,10,11,12,13,14,15]:
    Ec_part.append(Ec[i])

unknwns = [*x, *y[4:8]]

Sols = sp.linsolve(Ec_part,unknwns)

Sols = next( iter(Sols) )  

#sp.Matrix(sorted(Sols.items(), key=str))