I have made a solver which works fine for several other systems of equations that I've tried. However, my current system which has 21 equations (some non-linear and multivariate) when put into the same solver doesn't seem to work. The Newton's Method based matrix solver operates according to a while loop that compares the output residuals from guesses (called B2guessVariance) to a specified tolerance.
For some reason, and only with this example, I'm getting a strange output. After one loop iteration, my B2guessVariance matrix has a tau0 variable nested within the matrix. There is no variable in my code that uses the tau symbol. I suspect that it might have something to do with the gauss_jordan_solve operator I'm using from Sympy?
Some notes about the system in case it's bought up in any comments:
B2fsolveB2jsolveB2delta_x0B2guessVarianceThroughout the code, I have placed numerous "placeholders" to see where the code gets up to upon execution. I use sympy displays throughout my code but know that print(repr()) is preferable in most programs/terminals; for this reason the displays are commented out and the print(repr()) is used. The shortened code is as follows:
## 1.0 CODE SETTINGS, LIBRARIES & PHYSICAL CONSTANTS ##
from sympy.interactive import printing
printing.init_printing(use_latex = True)
import numpy as np
import sympy as sp
import math
import functools
from sympy import Matrix
from sympy.functions import sign
rho, g = sp.symbols('rho g')
cos = sp.cos
e = sp.exp
pi = np.pi
sqrt = sp.sqrt
## 2.0 PRV Matrix Constants ##
m_m, m_p, a_1, a_p, k_spr, theta = sp.symbols('m_m m_p a_1 a_p k_spr theta')
## 3.0 PRV FORCE BALANCE RELATIONSHIPS ##
xdot_m, a_2, q_3c, x_m, x_m0, t_0, t_f, t, P_sp = sp.symbols('xdot_m a_2 q_3c x_m x_m0 t_0 t_f t P_sp')
fa_2 = Matrix([a_2 - (1/(3700*(0.02732 - x_m)))]).doit()
fxdot_m = Matrix([0 - 3700*(0.02732 - x_m) * q_3c]).doit()
## 4.0 MAIN VALVE & PILOT VALVE FORCE BALANCES (m_m d^2m_m/dt^2) ##
h_in, h_out, h_c, q_m, x_p = sp.symbols('h_in h_out h_c q_m x_p')
fMValve = Matrix([rho * g * (h_in * a_1 + h_out * (a_2 - a_1) - h_c * a_2) - m_m * g * cos(theta) + (rho * ((q_m)**2)/a_1) - 0]).doit()
fPValve = Matrix([k_spr * (P_sp - x_p) - rho * g * h_out * a_p + m_p * g - 0]).doit()
## 5.0 FLOW EQUATIONS ##
C_vm, C_vfo, C_vp, C_vnvout, C_vnvin, h_t, q_1, q_2 = sp.symbols('C_vm C_vfo C_vp C_vnvout C_vnvin h_t q_1 q_2')
fC_vm = Matrix([0.02107 - 0.02962 * e(-51.1322 * x_m) + 0.0109 * e(-261 * x_m) - 0.00325 * e(-683.17 * x_m) + 0.0009 * e(-399.5 * x_m) - C_vm]).doit()
fq_m = Matrix([C_vm * (sqrt(((h_in - h_out)))) - q_m]).doit()
fq_1 = Matrix([C_vfo * (sqrt((h_in - h_t))) - q_1]).doit()
fC_vp = Matrix([0.0000753 * (1 - e(-1135 * x_p)) - C_vp]).doit()
fq_2 = Matrix([C_vp * (sqrt((h_t - h_out))) - q_2]).doit()
fmix = Matrix([q_1 + q_3c - q_2]).doit()
fq_3c = Matrix([h_c - h_t]).doit()
## 6.0 NETWORK EQUATIONS FROM PUMP TO ATMOSPHERIC OUTLET ##
C_valpha, C_vbeta, a_alpha, a_beta, h_pump, h_L1, h_L2, h_L3, h_L4, h_atm, f, L1, L2, L3, L4, D1, D2, D3, D4, A1, A2, A3, A4, v1, v2, v3, v4 = sp.symbols('C_valpha C_vbeta a_alpha a_beta h_pump h_L1 h_L2 h_L3 h_L4 h_atm f L1 L2 L3 L4 D1 D2 D3 D4 A1 A2 A3 A4 v1 v2 v3 v4')
#Pipe1
fp1q_m = Matrix([q_m - (A1 * v1)]).doit()
fh_L1 = Matrix([h_L1 - (f*(L1/D1)*((v1**2)/(2*g)))]).doit()
#AlphaValve
fCV1q_m = Matrix([C_valpha * a_alpha * (sqrt((h_pump - h_L1 - (h_in + h_L2)))) - q_m]).doit()
#Pipe2
fp2q_m = Matrix([q_m - (A2 * v2)]).doit()
fh_L2 = Matrix([h_L2 - (f*(L2/D2)*((v2**2)/(2*g)))]).doit()
#Pipe3
fp3q_m = Matrix([q_m - A3 * v3]).doit()
fh_L3 = Matrix([h_L3 - (f*(L3/D3)*((v3**2)/(2*g)))]).doit()
#BetaValve
fCV2q_m = Matrix([C_vbeta * a_beta * (sqrt((h_out - h_L3 - (h_atm + h_L4)))) - q_m]).doit()
#Pipe4
fp4q_m = Matrix([q_m - (A4 * v4)]).doit()
fh_L4 = Matrix([h_L4 - (f*(L4/D4)*((v4**2)/(2*g)))]).doit()
q_mscale, q_1scale, q_2scale, q_3cscale, x_mscale, x_pscale, P_spscale, a_2scale = sp.symbols('q_mscale q_1scale q_2scale q_3cscale x_mscale x_pscale P_spscale a_2scale')
## 7.1 INITIAL GUESSES & SOLVER CONDITIONS ##
guessrho, guessg, guessm_m, guessm_p, guessa_1, guessa_p, guessk_spr, guesstheta, guessP_sp, guessx_m, guessa_2, guessq_3c, guessh_in, guessh_out, guessh_c, guessq_m, guessx_p, guessC_vm, guessC_vfo, guessC_vp, guessh_t, guessq_1, guessq_2, guessC_valpha, guessC_vbeta, guessa_alpha, guessa_beta, guessh_pump, guessh_L1, guessh_L2, guessh_L3, guessh_L4, guessh_atm, guessf, guessL1, guessL2, guessL3, guessL4, guessD1, guessD2, guessD3, guessD4, guessA1, guessA2, guessA3, guessA4, guessv1, guessv2, guessv3, guessv4, guessq_mscale, guessq_1scale, guessq_2scale, guessq_3cscale, guessx_mscale, guessx_pscale, guessP_spscale, guessa_2scale = sp.symbols('guessrho guessg guessm_m guessm_p guessa_1 guessa_p guessk_spr guesstheta guessP_sp guessx_m guessa_2 guessq_3c guessh_in guessh_out guessh_c guessq_m guessx_p guessC_vm guessC_vfo guessC_vp guessh_t guessq_1 guessq_2 guessC_valpha guessC_vbeta guessa_alpha guessa_beta guessh_pump guessh_L1 guessh_L2 guessh_L3 guessh_L4 guessh_atm guessf guessL1 guessL2 guessL3 guessL4 guessD1 guessD2 guessD3 guessD4 guessA1 guessA2 guessA3 guessA4 guessv1 guessv2 guessv3 guessv4 guessq_mscale guessq_1scale guessq_2scale guessq_3cscale guessx_mscale guessx_pscale guessP_spscale guessa_2scale')
guessrho = 997 # physical constant
guessg = 9.81 # physical constant
guessm_m = 8 # physical constant
guessm_p = 0.1 # physical constant
guessa_1 = 0.0078 # physical constant
guessa_p = 0.00196 # physical constant
guessk_spr = 70000 # physical constant
guesstheta = 0 #0.000000000000000000000001 # physical constant 31/5/21 CHANGED SO IT'S NON-ZERO
guessP_sp = 0.00700000000
guessx_m = 0.013647064
guessa_2 = 0.01977
guessq_3c = -2.578E-24 #0.000000000000000000000001 # physical constant 31/5/21 CHANGED SO IT'S NON-ZERO
guessh_in = 41.38
guessh_out = 23.465726
guessh_c = 30.65
guessq_m = 0.02811
guessx_p = 0.0005878
guessC_vm = 0.006642
guessC_vfo = 0.00003
guessC_vp = 0.00003666
guessh_t = 30.65
guessq_1 = 0.00009826
guessq_2 = 0.00009826
guessC_valpha = 1
guessC_vbeta = 0.6
guessa_alpha = 0.01
guessa_beta = 0.01
guessh_pump = 50
guessh_L1 = 0.1959
guessh_L2 = 0.5224
guessh_L3 = 0.6529
guessh_L4 = 0.8619
guessh_atm = 0 #0.000000000000000000000001 # 31/5/21 CHANGED SO IT'S NON-ZERO
guessf = 0.02
guessL1 = 1.5
guessL2 = 4
guessL3 = 5
guessL4 = 6.6
guessD1 = 0.1
guessD2 = 0.1
guessD3 = 0.1
guessD4 = 0.1
guessA1 = 0.007854
guessA2 = 0.007854
guessA3 = 0.007854
guessA4 = 0.007854
guessv1 = 3.579
guessv2 = 3.579
guessv3 = 3.579
guessv4 = 3.579
guessq_mscale = 101.2 #q_m * 3600
guessq_1scale = 0.3537 #q_1 * 3600
guessq_2scale = 0.3537 #q_2 * 3600
guessq_3cscale = -9.281E-21 #q_3c * 36000
guessx_mscale = 13.6470637 #x_m * 1000
guessx_pscale = 0.5878 #x_p * 1000
guessP_spscale = 7 #P_sp * 1000
guessa_2scale = 197.7 #a_2 * (100^2)
#Broad Loop Solver Conditions
tolerance = 0.05
max_iter = 10
#!-!#! 7.2.2 BLOCK 2 (MULTIPLE UNKNOWN EQUATIONS) (THE PROBLEMS ARE WITH THIS BLOCK!!!!) #!-!#!-!#!-
print('*PROBLEM* SECTION 7.2.2 BLOCK 2 (MULTIPLE UNKNOWN EQUATIONS)')
B2FunctionMatrix = Matrix([[fxdot_m], [fa_2], [fMValve], [fPValve], [fq_3c], [fC_vm], [fq_m], [fq_1], [fC_vp], [fq_2], [fmix], [fp1q_m], [fh_L1], [fCV1q_m], [fp2q_m], [fh_L2], [fp3q_m], [fh_L3], [fCV1q_m], [fp4q_m], [fh_L4]]).subs(([rho, guessrho], [g, guessg], [m_m, guessm_m], [m_p, guessm_p], [a_1, guessa_1], [a_p, guessa_p], [k_spr, guessk_spr], [theta, guesstheta], [P_sp, guessP_sp], [C_vfo, guessC_vfo], [C_valpha, guessC_valpha], [C_vbeta, guessC_vbeta], [a_alpha, guessa_alpha], [a_beta, guessa_beta], [h_pump, guessh_pump], [h_atm, guessh_atm], [f, guessf], [L1, guessL1], [L2, guessL2], [L3, guessL3], [L4, guessL4], [D1, guessD1], [D2, guessD2], [D3, guessD3], [D4, guessD4], [A1, guessA1], [A2, guessA2], [A3, guessA3], [A4, guessA4], [P_spscale, guessP_spscale])).doit()
B2VariableMatrix = Matrix([[x_m], [a_2], [q_3c], [h_in], [h_out], [h_c], [q_m], [x_p], [C_vm], [C_vp], [h_t], [q_1], [q_2], [h_L1], [h_L2], [h_L3], [h_L4], [v1], [v2], [v3], [v4]]).doit()
#print('Block 2 (Multiple Unknown Equations)')
dof2 = abs(B2VariableMatrix.shape[0] - B2FunctionMatrix.shape[0])
if dof2 == 0:
B2guessmatrix = B2VariableMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq_3c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4])).doit()
#display(B2guessmatrix)
#display(B2FunctionMatrix)
B2iter_n = 0
B2tolerance = (Matrix.ones(B2FunctionMatrix.shape[0], 1) * tolerance).doit()
B2guessVariance = (Matrix.ones(B2FunctionMatrix.shape[0], 1) * 10).doit()
B2JacobianMatrix = (B2FunctionMatrix.jacobian(B2VariableMatrix)).doit()
#display(B2JacobianMatrix)
while (B2guessVariance[0,0] >= B2tolerance[0,0]) and (B2guessVariance[1,0] >= B2tolerance[1,0]) and (B2guessVariance[2,0] >= B2tolerance[2,0]) and (B2guessVariance[3,0] >= B2tolerance[3,0]) and (B2guessVariance[4,0] >= B2tolerance[4,0]) and (B2guessVariance[5,0] >= B2tolerance[5,0]) and (B2guessVariance[6,0] >= B2tolerance[6,0]) and (B2guessVariance[7,0] >= B2tolerance[7,0]) and (B2guessVariance[8,0] >= B2tolerance[8,0]) and (B2guessVariance[9,0] >= B2tolerance[9,0]) and (B2guessVariance[10,0] >= B2tolerance[10,0]) and (B2guessVariance[11,0] >= B2tolerance[11,0]) and (B2guessVariance[12,0] >= B2tolerance[12,0]) and (B2guessVariance[13,0] >= B2tolerance[13,0]) and (B2guessVariance[14,0] >= B2tolerance[14,0]) and (B2guessVariance[15,0] >= B2tolerance[15,0]) and (B2guessVariance[16,0] >= B2tolerance[16,0]) and (B2guessVariance[17,0] >= B2tolerance[17,0]) and (B2guessVariance[18,0] >= B2tolerance[18,0]) and (B2guessVariance[19,0] >= B2tolerance[19,0]) and (B2guessVariance[20,0] >= B2tolerance[20,0]):
#print('placeholder 1 (displays B2fsolve, the function matrix with guess values substituted)')
B2fsolve = B2FunctionMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq_3c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4])).doit()
#print(repr(B2fsolve))
#display(B2fsolve)
#print('placeholder 2 (displays B2jsolve, the jacobian matrix with guess values substituted)')
B2jsolve = B2JacobianMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq_3c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4])).doit()
#print(repr(B2jsolve))
#display(B2jsolve)
#print('placeholder 3')
B2delta_x0, fv = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
#The code does not run past this point, please help!
#print(repr(B2delta_x0))
#display(B2delta_x0)
#print('placeholder 4')
B2guessmatrix = B2VariableMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq_3c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4])).doit()
B2nextguessmatrix = B2delta_x0 + B2guessmatrix
guessx_m = B2nextguessmatrix[0]
guessa_2 = B2nextguessmatrix[1]
guessq_3c = B2nextguessmatrix[2]
guessh_in = B2nextguessmatrix[3]
guessh_out = B2nextguessmatrix[4]
guessh_c = B2nextguessmatrix[5]
guessq_m = B2nextguessmatrix[6]
guessx_p = B2nextguessmatrix[7]
guessC_vm = B2nextguessmatrix[8]
guessC_vp = B2nextguessmatrix[9]
guessh_t = B2nextguessmatrix[10]
guessq_1 = B2nextguessmatrix[11]
guessq_2 = B2nextguessmatrix[12]
guessh_L1 = B2nextguessmatrix[13]
guessh_L2 = B2nextguessmatrix[14]
guessh_L3 = B2nextguessmatrix[15]
guessh_L4 = B2nextguessmatrix[16]
guessv1 = B2nextguessmatrix[17]
guessv2 = B2nextguessmatrix[18]
guessv3 = B2nextguessmatrix[19]
guessv4 = B2nextguessmatrix[20]
B2guessVariance = abs(B2nextguessmatrix - B2guessmatrix)
B2guessmatrix = B2nextguessmatrix
#print('placeholder 5')
print('These are the two matrices which are compared in the while loop')
#display(B2guessVariance)
#display(B2tolerance)
print(repr(B2guessVariance))
print(repr(B2tolerance))
B2iter_n += 1
print(f'{B2iter_n} iterations have been completed so far. Moving onto the next iteration...')
if B2iter_n >= max_iter:
print('The maximum Number of iterations has been reached')
break
if (B2guessVariance[0,0] < B2tolerance[0,0]) and (B2guessVariance[1,0] < B2tolerance[1,0]) and (B2guessVariance[2,0] < B2tolerance[2,0]) and (B2guessVariance[3,0] < B2tolerance[3,0]) and (B2guessVariance[4,0] < B2tolerance[4,0]) and (B2guessVariance[5,0] < B2tolerance[5,0]) and (B2guessVariance[6,0] < B2tolerance[6,0]) and (B2guessVariance[7,0] < B2tolerance[7,0]) and (B2guessVariance[8,0] < B2tolerance[8,0]) and (B2guessVariance[9,0] < B2tolerance[9,0]) and (B2guessVariance[10,0] < B2tolerance[10,0]) and (B2guessVariance[11,0] < B2tolerance[11,0]) and (B2guessVariance[12,0] < B2tolerance[12,0]) and (B2guessVariance[13,0] < B2tolerance[13,0]) and (B2guessVariance[14,0] < B2tolerance[14,0]) and (B2guessVariance[15,0] < B2tolerance[15,0]) and (B2guessVariance[16,0] < B2tolerance[16,0]) and (B2guessVariance[17,0] < B2tolerance[17,0]) and (B2guessVariance[18,0] < B2tolerance[18,0]) and (B2guessVariance[19,0] < B2tolerance[19,0]) and (B2guessVariance[20,0] < B2tolerance[20,0]):
print('The solution Matrix is')
print(repr(B2nextguessmatrix))
#display(B2nextguessmatrix)
print(f'Which was achieved after {B2iter_n} iterations with a tolerance of {tolerance}.')
print(f'Displayed as integers, the solutions for variable {x_m} and {a_2} converge at {sp.Float(B2nextguessmatrix[0])} and {sp.Float(B2nextguessmatrix[1])} as floats respectively.')
else:
print('The Equation set has no solution or the initial guesses are too far from the solution.')
elif dof2 !=0:
print(f'This system has {B2FunctionMatrix.shape[0]} equations and {B2VariableMatrix.shape[0]} variables which represents a d.o.f value of {dof2} which ≠ 0. Therefore, the system cannot be solved.')
The output which shows the relational matrices used for iteration in the while loop is as follows:
*PROBLEM* SECTION 7.2.2 BLOCK 2 (MULTIPLE UNKNOWN EQUATIONS)
These are the two matrices which are compared in the while loop
Matrix([
[ Abs(0.00392183668209074*tau0 - 8.02180713197279e-7)],
[ Abs(0.00566975402880878*tau0 - 4.35392954981073e-6)],
[ Abs(7.3945310403193e-25*tau0 - 2.5781512492912e-24)],
[ Abs(4.81682791725627*tau0 - 0.000426049277400864)],
[ Abs(0.617506604286677*tau0 - 0.000126658199850027)],
[ Abs(4.04275751354693*tau0 - 0.00184684018329051)],
[ Abs(0.007854*tau0 - 5.33999999999951e-7)],
[Abs(0.000169107863923268*tau0 - 4.93212155862604e-8)],
[ Abs(0.00263411194406529*tau0 - 7.57176496023193e-7)],
[ Abs(7.41675230154456e-6*tau0 - 3.67497300027552e-9)],
[ Abs(4.04275751354693*tau0 - 0.00184684018329051)],
[ Abs(3.54463776960747e-6*tau0 - 3.53001248327642e-9)],
[ Abs(3.54463776960747e-6*tau0 - 3.53001248327642e-9)],
[ Abs(0.109449541284404*tau0 - 4.00458715596186e-5)],
[ Abs(0.291865443425076*tau0 - 0.000106788990825613)],
[ Abs(0.364831804281346*tau0 - 3.34862385320545e-5)],
[ Abs(0.481577981651376*tau0 - 0.000116201834862384)],
[ 1.0*Abs(tau0)],
[ 1.0*Abs(tau0)],
[ 1.0*Abs(tau0)],
[ Abs(tau0)]])
Matrix([
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05]])
1 iterations have been completed so far. Moving onto the next iteration...
This gives the error:
TypeError Traceback (most recent call last)
<ipython-input-1-e583c0741f11> in <module>
145 #display(B2JacobianMatrix)
146
--> 147 while (B2guessVariance[0,0] >= B2tolerance[0,0]) and (B2guessVariance[1,0] >= B2tolerance[1,0]) and (B2guessVariance[2,0] >= B2tolerance[2,0]) and (B2guessVariance[3,0] >= B2tolerance[3,0]) and (B2guessVariance[4,0] >= B2tolerance[4,0]) and (B2guessVariance[5,0] >= B2tolerance[5,0]) and (B2guessVariance[6,0] >= B2tolerance[6,0]) and (B2guessVariance[7,0] >= B2tolerance[7,0]) and (B2guessVariance[8,0] >= B2tolerance[8,0]) and (B2guessVariance[9,0] >= B2tolerance[9,0]) and (B2guessVariance[10,0] >= B2tolerance[10,0]) and (B2guessVariance[11,0] >= B2tolerance[11,0]) and (B2guessVariance[12,0] >= B2tolerance[12,0]) and (B2guessVariance[13,0] >= B2tolerance[13,0]) and (B2guessVariance[14,0] >= B2tolerance[14,0]) and (B2guessVariance[15,0] >= B2tolerance[15,0]) and (B2guessVariance[16,0] >= B2tolerance[16,0]) and (B2guessVariance[17,0] >= B2tolerance[17,0]) and (B2guessVariance[18,0] >= B2tolerance[18,0]) and (B2guessVariance[19,0] >= B2tolerance[19,0]) and (B2guessVariance[20,0] >= B2tolerance[20,0]):
148 #print('placeholder 1 (displays B2fsolve, the function matrix with guess values substituted)')
149 B2fsolve = B2FunctionMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq_3c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4])).doit()
C:\ProgramData\Anaconda3\lib\site-packages\sympy\core\relational.py in __bool__(self)
396
397 def __bool__(self):
--> 398 raise TypeError("cannot determine truth value of Relational")
399
400 def _eval_as_set(self):
TypeError: cannot determine the truth value of Relational
Obviously, as the tau0 symbol is in the output B2guessVariance matrix, it can't be related to the tolerance matrix for the while loop to continue iterating. Is there a way I can remove the tau0 variable so that the solver can execute properly and I can fix this?
After some research, this seems to be fixed by prespecifying the gauss_jordan_solve with the params return line and then replacing the taus by zeros.
So, this was resolved by changing from:
B2delta_x0, fv = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
to:
B2delta_x0, params = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
B2delta_x0 = B2delta_x0.xreplace({ tau:0 for tau in params })
Unfortunately, I don't really understand the theory behind as to why this works which is frustrating. But it's working for now and giving me good answers which is great!