CODEWITHSUNDEEP
pythonmatrixtypeerrorsympynewtons-method
Hendrix13
Hendrix13

Reputation: 125

Why do Sympy signum functions not work in Jacobian matrices when solving a system of equations?

I am trying to solve a system of 2 non-linear equations. These equations solve perfectly fine on their own, however, when I add Sympy's sign function to either of them they are unable to solve. The original system is as follows:

The modified system (identical to the original system but with a signum added to Equation 1) is:

The solver I am using is one I've created myself and it works fine for several other systems of equations that I've tried. The Newton's Method based matrix solver operates according to a while loop that compares the output residuals from guesses (called A0guessVariance) to a specified tolerance.

Some notes about the system in case it's bought up in any comments:

Here is the original minimised code without the signum function:

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.')

When running my code normally without the signum function the correct answer, (x_1, x_2) = (127, 54) is achieved as can be seen in the output below:

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.

However, when Equation 1 with the signum function is added the solver cannot execute. The final answer should cause the signum to return a positive value, I'm not sure why this is an issue. The same sample code but with the signum added to Equation 1 is below:

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.')

This returns the following output and error:

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

NOTE that the major difference between the first and second outputs seems to be A0JacobianMatrix at Placeholder 0. In the first output it is all numeric but in the second it is not; this seems like it could be a Sympy error. Any thoughts or help would be very much appreciated!

I have been getting a similar issue which I posted on Stack Overflow with a larger system of equations that has not been conclusively answered as of yet.

Upvotes: 0

Views: 724

Answers (1)

Hendrix13
Hendrix13

Reputation: 125

Big thanks to @OscarBenjamin for helping with the answer! (It's been a few days and I'm unable to mark his comment as the answer due to my low reputation I think?)

The issue was that the variables were not declared as real symbols. I presume Sympy assumes that the symbols are both real and complex by default.

To be clear, this issue is solved when every symbol is declared as real like this:

x_1, x_2 = sp.symbols('x_1 x_2', real = True)

Before it was simply like this:

x_1, x_2 = sp.symbols('x_1 x_2')

Upvotes: 0

Related Questions