问题描述
我有一个由 58 个方程组成的系统（其中一些方程是非线性的和/或多元的），我将它们分成多个块（通过 Trajan 算法），以便更容易地求解并加快计算速度。整个系统的解决方案肯定存在，因为我已经在软件（工程方程求解器 (EES)）上进行了尝试，并且得到了正确的答案。
现在，我正在尝试通过（主要是）Sympy 将我在 EES 上的原始系统重新创建/转换为 Python。现在，我接近完成，所有块都可以解决除了一个（这很令人抓狂！！！）。
有问题的代码与“Block 2”（B2）和牛顿法有关，其中：
- 方程/函数矩阵是“B2FunctionMatrix”，由 21 个方程组成。
- 变量、guess (initial)、next guess 和 Jacobian 矩阵在代码中都有相应的标签。
- 容差为 0.05
- 求解的函数矩阵为“B2fsolve”
- 求解的雅可比矩阵是“B2jsolve”
- 关于雅可比的求解函数矩阵是“B2delta_x0”
- 残差矩阵是“B2guessVariance”
在整个代码中，我放置了许多“占位符”以查看代码在执行时到达的位置。代码如下：
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],[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]))
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]])
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,guessq3_c],[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]))
B2iter_n = 0
B2tolerance = Matrix.ones(B2FunctionMatrix.shape[0],1) * tolerance
B2guessVariance = Matrix.ones(B2FunctionMatrix.shape[0],1) * 10
B2JacobianMatrix = B2FunctionMatrix.jacobian(B2VariableMatrix)
while B2guessVariance[0] > B2tolerance[0] or B2guessVariance[1] > B2tolerance[1]:
print('placeholder 1')
B2fsolve = B2FunctionMatrix.subs(([x_m,guessv4]))
print('placeholder 2')
B2jsolve = B2JacobianMatrix.subs(([x_m,guessv4]))
print('placeholder 3')
B2delta_x0,fv = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
print('placeholder 4')
B2guessmatrix = B2VariableMatrix.subs(([x_m,guessv4]))
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')
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] <= B2tolerance[0] or B2guessVariance[1] <= B2tolerance[1]:
print('The solution Matrix is')
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.')
我只看到代码中的“占位符 3”。需要明确的是，我得到的输出是：
Block 2 (Multiple UnkNown Equations)
placeholder 1
placeholder 2
placeholder 3
直到我收到以下让我非常困惑的错误：
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\cache.py in wrapper(*args,**kwargs)
71 try:
---> 72 retval = cfunc(*args,**kwargs)
73 except TypeError:
D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\operations.py in __new__(cls,evaluate,_sympify,*args)
84
---> 85 c_part,nc_part,order_symbols = cls.flatten(args)
86 is_commutative = not nc_part
D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\mul.py in flatten(cls,seq)
690 # order commutative part canonically
--> 691 _mulsort(c_part)
692
D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\mul.py in _mulsort(args)
30 # in-place sorting of args
---> 31 args.sort(key=_args_sortkey)
32
D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\basic.py in compare(self,other)
229 else:
--> 230 c = (l > r) - (l < r)
231 if c:
TypeError: unsupported operand type(s) for -: 'StrictGreaterThan' and 'StrictLessthan'
During handling of the above exception,another exception occurred:
TypeError Traceback (most recent call last)
<ipython-input-22-9ad901ca6789> in <module>
20 B2jsolve = B2JacobianMatrix.subs(([x_m,guessv4]))
21 print('placeholder 3')
---> 22 B2delta_x0,fv = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
23 print('placeholder 4')
24 B2guessmatrix = B2VariableMatrix.subs(([x_m,guessv4]))
D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\matrices.py in gauss_jordan_solve(self,B,freevar)
2158
2159 def gauss_jordan_solve(self,freevar=False):
-> 2160 return _gauss_jordan_solve(self,freevar=freevar)
2161
2162 def pinv_solve(self,arbitrary_matrix=None):
D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\solvers.py in _gauss_jordan_solve(M,freevar)
563
564 # solve by reduced row echelon form
--> 565 A,pivots = aug.rref(simplify=True)
566 A,v = A[:,:-B_cols],A[:,-B_cols:]
567 pivots = list(filter(lambda p: p < col,pivots))
D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\matrices.py in rref(self,iszerofunc,simplify,pivots,normalize_last)
169 normalize_last=True):
170 return _rref(self,iszerofunc=iszerofunc,simplify=simplify,--> 171 pivots=pivots,normalize_last=normalize_last)
172
173 echelon_form.__doc__ = _echelon_form.__doc__
D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\reductions.py in _rref(M,normalize_last)
304
305 mat,pivot_cols,_ = _row_reduce(M,simpfunc,--> 306 normalize_last,normalize=True,zero_above=True)
307
308 if pivots:
D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\reductions.py in _row_reduce(M,normalize_last,normalize,zero_above)
127 mat,swaps = _row_reduce_list(list(M),M.rows,M.cols,M.one,128 iszerofunc,normalize_last=normalize_last,--> 129 normalize=normalize,zero_above=zero_above)
130
131 return M._new(M.rows,mat),swaps
D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\reductions.py in _row_reduce_list(mat,rows,cols,one,zero_above)
67 pivot_offset,pivot_val,\
68 assumed_nonzero,newly_determined = _find_reasonable_pivot(
---> 69 get_col(piv_col)[piv_row:],simpfunc)
70
71 # _find_reasonable_pivot may have simplified some things
D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\determinant.py in _find_reasonable_pivot(col,simpfunc)
80 if possible_zeros[i] is not None:
81 continue
---> 82 simped = simpfunc(x)
83 is_zero = iszerofunc(simped)
84 if is_zero == True or is_zero == False:
D:\ProgramData\Anaconda3\lib\site-packages\sympy\simplify\simplify.py in simplify(expr,ratio,measure,rational,inverse,doit,**kwargs)
658 if expr.has(Piecewise):
659 # Fold into a single Piecewise
--> 660 expr = piecewise_fold(expr)
661 # Apply doit,if doit=True
662 expr = done(expr)
D:\ProgramData\Anaconda3\lib\site-packages\sympy\functions\elementary\piecewise.py in piecewise_fold(expr)
1127 args = expr.args
1128 # fold
-> 1129 folded = list(map(piecewise_fold,args))
1130 for ec in cartes(*[
1131 (i.args if isinstance(i,Piecewise) else
D:\ProgramData\Anaconda3\lib\site-packages\sympy\functions\elementary\piecewise.py in piecewise_fold(expr)
1132 [(i,true)]) for i in folded]):
1133 e,c = zip(*ec)
-> 1134 new_args.append((expr.func(*e),And(*c)))
1135
1136 return Piecewise(*new_args)
D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\cache.py in wrapper(*args,**kwargs)
72 retval = cfunc(*args,**kwargs)
73 except TypeError:
---> 74 retval = func(*args,**kwargs)
75 return retval
76
D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\operations.py in __new__(cls,*args)
83 return args[0]
84
---> 85 c_part,order_symbols = cls.flatten(args)
86 is_commutative = not nc_part
87 obj = cls._from_args(c_part + nc_part,is_commutative)
D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\mul.py in flatten(cls,seq)
689
690 # order commutative part canonically
--> 691 _mulsort(c_part)
692
693 # current code expects coeff to be always in slot-0
D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\mul.py in _mulsort(args)
29 def _mulsort(args):
30 # in-place sorting of args
---> 31 args.sort(key=_args_sortkey)
32
33
D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\basic.py in compare(self,other)
228 c = l.compare(r)
229 else:
--> 230 c = (l > r) - (l < r)
231 if c:
232 return c
TypeError: unsupported operand type(s) for -: 'StrictGreaterThan' and 'StrictLessthan'
其他 3 个块完全解决。我在想这可能与块中稀疏雅可比矩阵的大尺寸有关（“B2Jsolve”是一个稀疏的 21 x 21 方阵）？任何想法将不胜感激！感谢阅读:)
B2Jsolve 的输出如下： B2jsolve (This is a merged jpg file only with good resolution when you zoom in to make it clearer!)
编辑（添加源代码）：这是我添加了##headings 的代码的简化版本，以便更清晰一些。求解块位于代码的底部。我认为问题可能与对一些使用 Sympy 符号函数的方程（第 33、34、36、47 和 55 行的原始方程）取雅可比 (B2jsolve) 有关。
## 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)))])
fxdot_m = Matrix([0 - 3700*(0.02732 - x_m) * q_3c])
## 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])
fPValve = Matrix([k_spr * (P_sp - x_p) - rho * g * h_out * a_p + m_p * g - 0])
## 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])
fq_m = Matrix([C_vm * (sqrt((abs(h_in - h_out)))) * sign(h_in - h_out) - q_m])
fq_1 = Matrix([C_vfo * (sqrt(abs(h_in - h_t))) * sign(h_in - h_t) - q_1])
fC_vp = Matrix([0.0000753 * (1 - e(-1135 * x_p)) - C_vp])
fq_2 = Matrix([C_vp * (sqrt(abs(h_t - h_out))) * sign(h_t - h_out) - q_2])
fmix = Matrix([q_1 + q_3c - q_2])
fq_3c = Matrix([h_c - h_t])
## 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)])
fh_L1 = Matrix([h_L1 - (f*(L1/D1)*((v1**2)/(2*g)))])
#AlphaValve
fCV1q_m = Matrix([C_valpha * a_alpha * (sqrt(abs(h_pump - h_L1 - (h_in + h_L2)))) * sign(h_pump - h_L1 - (h_in + h_L2)) - q_m])
#Pipe2
fp2q_m = Matrix([q_m - (A2 * v2)])
fh_L2 = Matrix([h_L2 - (f*(L2/D2)*((v2**2)/(2*g)))])
#Pipe3
fp3q_m = Matrix([q_m - A3 * v3])
fh_L3 = Matrix([h_L3 - (f*(L3/D3)*((v3**2)/(2*g)))])
#BetaValve
fCV2q_m = Matrix([C_vbeta * a_beta * (sqrt(abs(h_out - h_L3 - (h_atm + h_L4)))) * sign(h_out - h_L3 - (h_atm + h_L4)) - q_m])
#Pipe4
fp4q_m = Matrix([q_m - (A4 * v4)])
fh_L4 = Matrix([h_L4 - (f*(L4/D4)*((v4**2)/(2*g)))])
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],[v4]])
print('Block 2 (Multiple UnkNown Equations)')
dof2 = abs(B2VariableMatrix.shape[0] - B2FunctionMatrix.shape[0])
if dof2 == 0:
B2guessmatrix = B2VariableMatrix.subs(([x_m,guessq_3c],guessv4]))
#display(B2guessmatrix)
#display(B2FunctionMatrix)
B2iter_n = 0
B2tolerance = Matrix.ones(B2FunctionMatrix.shape[0],1) * 10
B2JacobianMatrix = B2FunctionMatrix.jacobian(B2VariableMatrix)
#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,guessv4]))
print(repr(B2fsolve))
#display(B2fsolve)
print('placeholder 2 (displays B2jsolve,the jacobian matrix with guess values substituted)')
B2jsolve = B2JacobianMatrix.subs(([x_m,guessv4]))
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,guessv4]))
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')
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] < B2tolerance[1,0] < B2tolerance[2,0] < B2tolerance[3,0] < B2tolerance[4,0] < B2tolerance[5,0] < B2tolerance[6,0] < B2tolerance[7,0] < B2tolerance[8,0] < B2tolerance[9,0] < B2tolerance[10,0] < B2tolerance[11,0] < B2tolerance[12,0] < B2tolerance[13,0] < B2tolerance[14,0] < B2tolerance[15,0] < B2tolerance[16,0] < B2tolerance[17,0] < B2tolerance[18,0] < B2tolerance[19,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 system cannot be solved.')
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）