使用Newton-Raphson法求根时出现‘cannot determine truth value of Relational’错误的技术求助

问题描述 #

我在使用牛顿-拉夫逊法求解根时遇到了错误：cannot determine truth value of Relational。以下是我的代码，主要是用SymPy定义了一些符号表达式，然后尝试用牛顿法求解行列式的根：

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)

错误原因分析 #

这个错误的核心原因是你在牛顿-拉夫逊的循环中直接使用了SymPy的符号表达式进行数值比较，而不是先将符号表达式转换为数值结果：

1. 你定义的det和ddet是SymPy的符号表达式，依赖于符号变量x。但你的f(x)和g(x)函数只是直接返回这些符号表达式，没有代入输入的数值x进行计算。
2. 当你在g(x0) == 0.0或abs(f(x1)) > e中进行比较时，得到的是SymPy的Relational对象（比如Eq或Gt），Python无法直接判断这些符号关系的真假，因此抛出错误。
3. 另外，你的循环变量k和你定义的符号变量k重名了（symbols('x c1 c2 c3 k')），这会在sum(U[k] for k in range(N+1))中覆盖符号k，导致后续符号表达式出现问题。

解决方案 #

我们需要将符号表达式转换为可以计算数值的函数，并修正变量名冲突问题，具体步骤如下：

1. 修正变量名冲突

将符号变量k重命名为其他名字，比如k_sym，避免和循环变量k冲突：

# 原代码
# x,c1, c2, c3, k = symbols('x c1 c2 c3 k')
# 修改为
x,c1, c2, c3, k_sym = symbols('x c1 c2 c3 k_sym')

2. 将符号表达式转换为数值函数

使用sympy.lambdify将det和ddet转换为可以接收数值x并返回数值结果的函数：

# 替换原有的f和g函数定义
# 把符号表达式转换为数值函数，x是输入变量
f_num = lambdify(x, det, 'numpy')
g_num = lambdify(x, ddet, 'numpy')

# 定义用于牛顿法的f和g，接收数值x并返回数值
def f(x_val):
    return f_num(x_val)

def g(x_val):
    return g_num(x_val)

3. 确保循环中的条件是数值比较

在牛顿-拉夫逊的循环中，所有比较操作都基于数值结果，此时f(x1)和g(x0)都是数值，不会再出现符号关系的判断问题。

修正后的完整代码 #

import sympy as sp
import numpy as np
from scipy import io, integrate, linalg, signal
from scipy.sparse.linalg import cg, eigs

# 避免重复导入，删除重复的import语句
# 定义符号变量，将k改为k_sym避免和循环变量冲突
x,c1, c2, c3, k_sym = symbols('x c1 c2 c3 k_sym')
A11=4.9091
B11=0.0294
D11=0.0014
I0=140.8182
I1=0.1191
I2=0.0307

# 定义e1到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)

# 初始化U和W列表
U = [0 for _ in range(50)]
W = [0 for _ in range(50)]

# 边界条件
U[0] = 0
U[1] = c1
W[0] = 0
W[1] = 0
W[2] = c2
W[3] = c3

# 循环计算U和W的项
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))

    # 计算求和方程，注意循环变量用k，和符号k_sym不冲突
    f1 = sum(U[k] for k in range(N+1))
    f2 = sum(W[k] for k in range(N+1))
    f3 = sum(k * W[k] for k in range(N+1))
    
    f11 = sp.expand(f1)
    f22 = sp.expand(f2)
    f33 = sp.expand(f3)

    # 整理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]

    # 构建系数矩阵A和常数向量b
    A = sp.Matrix([[eqt.coeff(c1) for eqt in system],
                   [eqt.coeff(c2) for eqt in system],
                   [eqt.coeff(c3) for eqt in system]])

    b = sp.Matrix([eqt.subs({c1: 0, c2: 0, c3: 0}) for eqt in system])

    # 计算行列式和导数
    det = A.det()
    ddet = sp.diff(det, x)
    
    # 将符号表达式转换为数值函数
    f_num = sp.lambdify(x, det, 'numpy')
    g_num = sp.lambdify(x, ddet, 'numpy')

    def f(x_val):
        return f_num(x_val)

    def g(x_val):
        return g_num(x_val)
    
    # 牛顿-拉夫逊实现
    def newtonRaphson(x0, e, Nst):
        print('\n\n*** NEWTON RAPHSON METHOD IMPLEMENTATION ***')
        step = 1
        x_current = float(x0)
        while step <= Nst:
            g_val = g(x_current)
            if abs(g_val) < 1e-10:  # 避免除以0，用极小值判断
                print('Divide by zero error! Derivative is too small.')
                break
            
            f_val = f(x_current)
            x_next = x_current - f_val / g_val
            f_next = f(x_next)
            
            # 判断是否收敛
            if abs(f_next) <= e:
                print(f'\nRequired root is: {x_next:.8f}')
                break
            
            x_current = x_next
            step += 1
        
        if step > Nst:
            print('\nNot Convergent.')


    # 输入参数
    x0 = 0.1
    e = 0.0001
    Nst = 10

    # 启动牛顿法
    newtonRaphson(x0, e, Nst)

额外说明 #

• lambdify函数可以将SymPy符号表达式转换为NumPy兼容的数值函数，这样计算效率更高，也避免了符号操作的问题。
• 我还调整了牛顿法的循环逻辑，让它更简洁，同时添加了对导数极小值的判断，避免除以0的情况。
• 变量名冲突是很容易被忽略的问题，一定要确保循环变量和符号变量不重名，否则会破坏符号表达式的结构。

备注：内容来源于stack exchange，提问作者YouTldi