python - 为什么在求解方程组时 Sympy 符号函数在雅可比矩阵中不起作用？
问题描述
我正在尝试求解一个由 2 个非线性方程组成的系统。这些方程本身就可以很好地求解，但是，当我将 Sympy 的sign函数添加到它们中的任何一个时，它们都无法求解。原系统如下：
- 等式1；x_2 = 0.00005x_1^3 + 0.0035x_1^2 + 0.6134x_1 - 183.6
- 等式2；x_2 = -0.0005x_1^2 - 0.0028x_1 + 62.765
修改后的系统（与原始系统相同，但在等式 1 中添加了一个符号）是：
- 等式1；x_2 = (0.00005x_1^3 + 0.0035x_1^2 + 0.6134x_1 - 183.6) * 符号(x_1 - x_2)
- 等式2；x_2 = -0.0005x_1^2 - 0.0028x_1 + 62.765
我使用的求解器是我自己创建的，它适用于我尝试过的其他几个方程组。基于牛顿法的矩阵求解器根据将猜测的输出残差（称为A0guessVariance）与指定容差进行比较的 while 循环进行操作。
有关该系统的一些说明，以防在任何评论中被购买：
- 我的 2 方程系统的解决方案肯定存在
- 求解的函数矩阵为
A0fsolve - 求解的雅可比矩阵是
A0JacobianMatrix。用猜测值代替它是A0jsolve - 关于雅可比的求解函数矩阵是
A0delta_x0 - 残差矩阵是
A0guessVariance - 无论输入的方程是线性的、非线性的和/或多元的等，我制作的求解器以前都可以工作。
这是没有符号函数的原始最小化代码：
import sympy as sp
from sympy.interactive import printing
printing.init_printing(use_latex = True)
from sympy import Matrix
from sympy.functions import sign
x_1, x_2 = sp.symbols('x_1 x_2')
#Input each equation into the "A0Function Matrix"
Equation1 = Matrix([((5*10**-5)*x_1**3 + +0.0035*x_1**2 + 0.6134*x_1 - 183.6 - x_2)])
#Equation1 = Matrix([((5*10**-5)*x_1**3 + +0.0035*x_1**2 + 0.6134*x_1 - 183.6 - x_2) * sign(x_1 - x_2)]) #UNCOMMENT FOR ERROR
Equation2 = Matrix([-0.0005*x_1**2 - 0.0028*x_1 + 62.765 - x_2])
A0FunctionMatrix = Matrix([Equation1, Equation2])
print('The System Function Matrix is')
#display(A0FunctionMatrix)
print(repr(A0FunctionMatrix))
#Variable Definitions
A0VariableMatrix = Matrix([x_1, x_2])
# The initial guesses must first be defined
guessx_1, guessx_2 = sp.symbols('guessx_1 guessx_2')
guessx_1 = 125
guessx_2 = 50
#Broad Loop Solver Conditions
tolerance = 0.05
max_iter = 5
print('SOLVER:')
dof = abs(A0VariableMatrix.shape[0] - A0FunctionMatrix.shape[0])
if dof == 0:
A0guessmatrix = A0VariableMatrix.subs(([x_1, guessx_1], [x_2, guessx_2]))
A0iter_n = 0
A0tolerance = Matrix.ones(A0FunctionMatrix.shape[0], 1) * tolerance
A0guessVariance = Matrix.ones(A0FunctionMatrix.shape[0], 1) * 10
A0JacobianMatrix = A0FunctionMatrix.jacobian(A0VariableMatrix)
print('placeholder 0 (displays the original Jacobian Matrix)')
print(repr(A0JacobianMatrix))
while (A0guessVariance[0,0] > A0tolerance[0,0]) and (A0guessVariance[1,0] > A0tolerance[1,0]):
#print('Placeholder 1 (displays "A0fsolve", the function matrix with guess values substituted)')
A0fsolve = A0FunctionMatrix.subs(([x_1, guessx_1], [x_2, guessx_2]))
#display(A0fsolve)
#print(repr(A0fsolve))
print('Placeholder 2 (displays "A0jsolve", the jacobian matrix with guess values substituted)')
A0jsolve = A0JacobianMatrix.subs(([x_1, guessx_1], [x_2, guessx_2]))
#display(A0jsolve)
print(repr(A0jsolve))
A0delta_x0, fv = A0jsolve.gauss_jordan_solve(-1*A0fsolve)
A0guessmatrix = A0VariableMatrix.subs(([x_1, guessx_1], [x_2, guessx_2]))
A0nextguessmatrix = A0delta_x0 + A0guessmatrix
guessx_1 = A0nextguessmatrix[0]
guessx_2 = A0nextguessmatrix[1]
A0guessVariance = abs(A0nextguessmatrix - A0guessmatrix)
A0guessmatrix = A0nextguessmatrix
A0iter_n += 1
print(f'{A0iter_n} iterations have been completed so far. Moving onto the next iteration...')
if A0iter_n >= max_iter:
print('The maximum Number of iterations has been reached')
break
if (A0guessVariance[0,0] <= A0tolerance[0,0]) and (A0guessVariance[1,0] <= A0tolerance[1,0]):
print('The solution Matrix is')
#display(A0nextguessmatrix)
print(repr(A0nextguessmatrix))
print(f'Which was achieved after {A0iter_n} iterations with a tolerance of {tolerance}.')
print(f'Displayed as integers, the solutions for variable x and y converge at {sp.Float(A0nextguessmatrix[0])} and {sp.Float(A0nextguessmatrix[1])} as floats respectively.')
else:
print('The Equation set has no solution or the initial guesses are too far from the solution.')
elif dof !=0:
print(f'This system has {A0FunctionMatrix.shape[0]} equations and {A0VariableMatrix.shape[0]} variables which represents a d.o.f value of {dof} which ≠ 0. Therefore, the system cannot be solved.')
在没有符号函数的情况下正常运行我的代码时，正确答案是 (x_1, x_2) = (127, 54) ，如下面的输出所示：
The System Function Matrix is
Matrix([
[5.0e-5*x_1**3 + 0.0035*x_1**2 + 0.6134*x_1 - x_2 - 183.6],
[ -0.0005*x_1**2 - 0.0028*x_1 - x_2 + 62.765]])
SOLVER:
placeholder 0 (displays the original Jacobian Matrix)
Matrix([
[0.00015*x_1**2 + 0.007*x_1 + 0.6134, -1],
[ -0.001*x_1 - 0.0028, -1]])
Placeholder 2 (displays "A0jsolve", the jacobian matrix with guess values substituted)
Matrix([
[3.83215, -1],
[-0.1278, -1]])
1 iterations have been completed so far. Moving onto the next iteration...
Placeholder 2 (displays "A0jsolve", the jacobian matrix with guess values substituted)
Matrix([
[ 3.93615930824608, -1],
[-0.130119158070178, -1]])
2 iterations have been completed so far. Moving onto the next iteration...
The solution Matrix is
Matrix([
[127.288913112303],
[54.3073578000086]])
Which was achieved after 2 iterations with a tolerance of 0.05.
Displayed as integers, the solutions for variable x and y converge at 127.288913112303 and 54.3073578000086 as floats respectively.
但是，当添加带有符号函数的公式 1时，求解器将无法执行。最终答案应该导致符号返回正值，我不确定为什么这是一个问题。下面是相同的示例代码，但在公式 1 中添加了符号：
import sympy as sp
from sympy.interactive import printing
printing.init_printing(use_latex = True)
from sympy import Matrix
from sympy.functions import sign
x_1, x_2 = sp.symbols('x_1 x_2')
#Input each equation into the "A0Function Matrix"
#Equation1 = Matrix([((5*10**-5)*x_1**3 + +0.0035*x_1**2 + 0.6134*x_1 - 183.6 - x_2)])
Equation1 = Matrix([((5*10**-5)*x_1**3 + +0.0035*x_1**2 + 0.6134*x_1 - 183.6 - x_2) * sign(x_1 - x_2)]) #UNCOMMENT FOR ERROR
Equation2 = Matrix([-0.0005*x_1**2 - 0.0028*x_1 + 62.765 - x_2])
A0FunctionMatrix = Matrix([Equation1, Equation2])
print('The System Function Matrix is')
#display(A0FunctionMatrix)
print(repr(A0FunctionMatrix))
#Variable Definitions
A0VariableMatrix = Matrix([x_1, x_2])
# The initial guesses must first be defined
guessx_1, guessx_2 = sp.symbols('guessx_1 guessx_2')
guessx_1 = 125
guessx_2 = 50
#Broad Loop Solver Conditions
tolerance = 0.05
max_iter = 5
print('SOLVER:')
dof = abs(A0VariableMatrix.shape[0] - A0FunctionMatrix.shape[0])
if dof == 0:
A0guessmatrix = A0VariableMatrix.subs(([x_1, guessx_1], [x_2, guessx_2]))
A0iter_n = 0
A0tolerance = Matrix.ones(A0FunctionMatrix.shape[0], 1) * tolerance
A0guessVariance = Matrix.ones(A0FunctionMatrix.shape[0], 1) * 10
A0JacobianMatrix = A0FunctionMatrix.jacobian(A0VariableMatrix)
print('placeholder 0 (displays the original Jacobian Matrix)')
print(repr(A0JacobianMatrix))
while (A0guessVariance[0,0] > A0tolerance[0,0]) and (A0guessVariance[1,0] > A0tolerance[1,0]):
#print('Placeholder 1 (displays "A0fsolve", the function matrix with guess values substituted)')
A0fsolve = A0FunctionMatrix.subs(([x_1, guessx_1], [x_2, guessx_2]))
#display(A0fsolve)
#print(repr(A0fsolve))
print('Placeholder 2 (displays "A0jsolve", the jacobian matrix with guess values substituted)')
A0jsolve = A0JacobianMatrix.subs(([x_1, guessx_1], [x_2, guessx_2]))
#display(A0jsolve)
print(repr(A0jsolve))
A0delta_x0, fv = A0jsolve.gauss_jordan_solve(-1*A0fsolve)
A0guessmatrix = A0VariableMatrix.subs(([x_1, guessx_1], [x_2, guessx_2]))
A0nextguessmatrix = A0delta_x0 + A0guessmatrix
guessx_1 = A0nextguessmatrix[0]
guessx_2 = A0nextguessmatrix[1]
A0guessVariance = abs(A0nextguessmatrix - A0guessmatrix)
A0guessmatrix = A0nextguessmatrix
A0iter_n += 1
print(f'{A0iter_n} iterations have been completed so far. Moving onto the next iteration...')
if A0iter_n >= max_iter:
print('The maximum Number of iterations has been reached')
break
if (A0guessVariance[0,0] <= A0tolerance[0,0]) and (A0guessVariance[1,0] <= A0tolerance[1,0]):
print('The solution Matrix is')
#display(A0nextguessmatrix)
print(repr(A0nextguessmatrix))
print(f'Which was achieved after {A0iter_n} iterations with a tolerance of {tolerance}.')
print(f'Displayed as integers, the solutions for variable x and y converge at {sp.Float(A0nextguessmatrix[0])} and {sp.Float(A0nextguessmatrix[1])} as floats respectively.')
else:
print('The Equation set has no solution or the initial guesses are too far from the solution.')
elif dof !=0:
print(f'This system has {A0FunctionMatrix.shape[0]} equations and {A0VariableMatrix.shape[0]} variables which represents a d.o.f value of {dof} which ≠ 0. Therefore, the system cannot be solved.')
这将返回以下输出和错误：
The System Function Matrix is
Matrix([
[(5.0e-5*x_1**3 + 0.0035*x_1**2 + 0.6134*x_1 - x_2 - 183.6)*sign(x_1 - x_2)],
[ -0.0005*x_1**2 - 0.0028*x_1 - x_2 + 62.765]])
SOLVER:
placeholder 0 (displays the original Jacobian Matrix)
Matrix([
[(0.00015*x_1**2 + 0.007*x_1 + 0.6134)*sign(x_1 - x_2) + (5.0e-5*x_1**3 + 0.0035*x_1**2 + 0.6134*x_1 - x_2 - 183.6)*Derivative(sign(x_1 - x_2), x_1), (5.0e-5*x_1**3 + 0.0035*x_1**2 + 0.6134*x_1 - x_2 - 183.6)*Derivative(sign(x_1 - x_2), x_2) - sign(x_1 - x_2)],
[ -0.001*x_1 - 0.0028, -1]])
Placeholder 2 (displays "A0jsolve", the jacobian matrix with guess values substituted)
Matrix([
[3.83215 - 4.58125*Subs(Derivative(sign(x_1 - 50), x_1), x_1, 125), -4.58125*Subs(Derivative(sign(125 - x_2), x_2), x_2, 50) - 1],
[ -0.1278, -1]])
1 iterations have been completed so far. Moving onto the next iteration...
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-13-193820857f4b> in <module>
39 print(repr(A0JacobianMatrix))
40
---> 41 while (A0guessVariance[0,0] > A0tolerance[0,0]) and (A0guessVariance[1,0] > A0tolerance[1,0]):
42 #print('Placeholder 1 (displays "A0fsolve", the function matrix with guess values substituted)')
43 A0fsolve = A0FunctionMatrix.subs(([x_1, guessx_1], [x_2, guessx_2]))
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 truth value of Relational
请注意，第一个和第二个输出之间的主要区别似乎A0JacobianMatrix在占位符 0 处。在第一个输出中，它都是数字，但在第二个输出中不是；这似乎是一个 Sympy 错误。任何想法或帮助将不胜感激！
我收到了一个类似的问题，我在 Stack Overflow 上发布了一个更大的方程组，到目前为止还没有最终答案。
解决方案
非常感谢@OscarBenjamin 帮助解答！（已经有几天了，由于我认为我的声誉很低，我无法将他的评论标记为答案？）
问题是变量没有被声明为真正的符号。我假设 Sympy 默认假设符号既真实又复杂。
需要明确的是，当每个符号都像这样声明为真实时，这个问题就解决了：
x_1, x_2 = sp.symbols('x_1 x_2', real = True)
之前简直就是这样：
x_1, x_2 = sp.symbols('x_1 x_2')
推荐阅读
- javascript - 使用 MongoDB 和 Nodejs 创建嵌套的子文档/子元素
- xml - 如何使用 PowerShell 从 XML 文件中获取退出代码？
- swift - 如何在 Xcode swift 中以编程方式更改 UIBarButtunItem 功能
- postfix-mta - 如何将 Postfix 配置为仅中继来自特定域的电子邮件？
- python-3.x - 不止一种方法对硒元素进行定位
- istio - 带有 Keycloak 的 Istio OAuth2
- java - Springboot依赖注入在休眠实体中失败
- android - FireBase 卡在 Android 的“运行您的应用程序以验证安装”上
- java - 如何使用 Sikuli 终止应用程序
- php - Laravel 在 x 小时分钟秒后延迟队列作业调度