I'm using the Newton-Raphson method for root-finding and ran into the error: 'cannot determine truth value of Relational.' I'd really appreciate any insights into what might be causing this. Thanks!
import sympy as sp
import numpy as np
from scipy import io, integrate, linalg, signal
from scipy import io, integrate, linalg, signal
from scipy.sparse.linalg import cg, eigs
import sympy as sp
import math
from math import *
from sympy import *
from sympy import symbols
from sympy import lambdify
# Define symbolic variable
# Define symbols for c1, c2, c3, and other variables
x,c1, c2, c3, k = symbols('x c1 c2 c3 k')
A11=4.9091
B11=0.0294
D11=0.0014
I0=140.8182
I1=0.1191
I2=0.0307
# Define expressions for e1, e2, e3, e4, e5, and e6
e1 = -I0 * x / A11
e2 = I1 * x / A11
e3 = B11 / A11
e4 = -(B11 * e1 + I1 * x) / (B11 * e3 - D11)
e5 = -I0 * x / (B11 * e3 - D11)
e6 = (-B11 * e2 + I2 * x) / (B11 * e3 - D11)
# Define U and W as symbolic lists
U = [0 for _ in range(50)] # Extend this as needed
W = [0 for _ in range(50)] # Extend this as needed
# Initialize the given boundary conditions
U[0] = 0
U[1] = c1
W[0] = 0
W[1] = 0
W[2] = c2
W[3] = c3
# Loop over N and K, as in the Maple code
for N in range(12, 30, 2):
for K in range(0, N):
U[K+2] = (e1*U[K] + e2*(K+1)*W[K+1] + e3*(K+3)*(K+2)*(K+1)*W[K+3]) / ((K+2)*(K+1))
W[K+4] = (e4*(K+1)*U[K+1] + e5*W[K] + e6*(K+2)*(K+1)*W[K+2]) / ((K+4)*(K+3)*(K+2)*(K+1))
# Define the sum equations
f1 = sum(U[k] for k in range(N+1)) # Sum for U
f2 = sum(W[k] for k in range(N+1)) # Sum for W
f3 = sum(k * W[k] for k in range(N+1)) # Sum for W weighted by k
f11=sp.expand(f1)
f22=sp.expand(f2)
f33=sp.expand(f3)
# Collect terms involving c1, c2, c3
eqt1 = sp.collect(f11, {c1, c2, c3})
eqt2 = sp.collect(f22, {c1, c2, c3})
eqt3 = sp.collect(f33, {c1, c2, c3})
system = [eqt1, eqt2, eqt3]
# Extract the coefficients of c1, c2, c3
A = sp.Matrix([[eqt.coeff(c1) for eqt in system],
[eqt.coeff(c2) for eqt in system],
[eqt.coeff(c3) for eqt in system]])
# Create the vector for the constants
b = sp.Matrix([eqt.subs({c1: 0, c2: 0, c3: 0}) for eqt in system])
# Compute the determinant of the matrix A
det = A.det()
ddet=diff(det,x)
# Defining Function
def f(x):
return det
# Defining derivative of function
def g(x):
return ddet
# Implementing Newton Raphson Method
def newtonRaphson(x0,e,Nst):
print('\n\n*** NEWTON RAPHSON METHOD IMPLEMENTATION ***')
step = 1
flag = 1
condition = True
while condition:
if g(x0) == 0.0:
print('Divide by zero error!')
break
x1 = x0 - f(x0)/g(x0)
#print('Iteration-%d, x1 = %0.6f and f(x1) = %0.6f' % (step, x1, f(x1)))
x0 = x1
step = step + 1
if step > Nst:
flag = 0
break
condition = abs(f(x1)) > e
if flag==1:
print('\nRequired root is: %0.8f' % x1)
else:
print('\nNot Convergent.')
# Input Section
x0 =0.1 #input('Enter Guess: ')
e = 0.0001#input('Tolerable Error: ')
Nst =10# input('Maximum Step: ')
# Converting x0 and e to float
x0 = float(x0)
e = float(e)
# Converting N to integer
Nst = int(Nst)
#Note: You can combine above three section like this
# x0 = float(input('Enter Guess: '))
# e = float(input('Tolerable Error: '))
# N = int(input('Maximum Step: '))
# Starting Newton Raphson Method
newtonRaphson(x0,e,Nst)
I'm using the Newton-Raphson method for root-finding and ran into the error: 'cannot determine truth value of Relational.'
There is a deep misunderstanding about what sympy symbols are.
sympy.symbols('x') has nothing to do with python variable x. Except that, most of the times, we x=sympy.symbols('x'), so that python variable x is a sympy object representing symbol x.
Which is what you did.
After
x,c1, c2, c3, k = symbols('x c1 c2 c3 k')
variable x does contain the symbol x
But inside code of f
def f(x):
return det
variable x is the argument of function f.
While symbol x (that is used by symbolic expression det) has nothing to do with that.
Let's use a simpler example.
What you do is the equivalent of
x=sympy.symbols('x')
det=2*x+1
def f(x):
return det
print(f(12))
Obviously expecting that it would return 25. But it won't.
det is symbolic expression 2x+1. So f returns symbolic expression 2x+1 regardless of whatever argument you passed to it, that is simply never used here. The fact that you choose to name parameter of f x is irrelevant here.
What you probably want to do is
x=sympy.symbols('x')
det=2*x+1
f = sympy.lambdify(x, det)
print(f(12))
lambdify will build a function whose parameter is the passed symbol (or symbols if you pass more than one), and return value is the evaluation of the expression (det) when that symbol is evaluated to the values passed as argument of the said built function.
So, in short, here f(12) will indeed return 25.
Note again that here, it is the symbol that matters, no the name of the variable that contains it.
Let's rephrase that example
varX=sympy.symbols('symx')
det=2*varX+1
# det is 2symx+1
f = sympy.lambdify(varX, det)
# f is a function f(argx) that returns the value of expression det if
# symx is replaced by argx
print(f(12))
# 25
So, in your case
f=lambdify(x, det)
g=lambdify(x, ddet)
Not claiming that it will work (there seems to be other numerical problems). But it solve your current problem
(I note that you imported lambdify, but never used it. So I guess that what I say here should not be entirely new to you)