Python牛顿-拉夫逊迭代代码出现除零错误求助

问题概述 #

编写用于数值计算的牛顿-拉夫逊迭代Python程序时，触发division by zero错误，错误出现在f(x)函数的返回行。尝试添加x == 0的判断逻辑后，问题仍未解决。

完整原始代码 #

# Newton-Raphson iterative scheme
qin_secs = 100 
W = 3
L = 7000
A = W * L
n = 0.04
S = 0.007
b = 100
deltatime = 1

# define the function f(x)
def f(x):
    return qin_secs-(((x-b)/deltatime)*A)-((W*x)/n)*((W*x)/((W+2)*x))**0.66*(S)**0.5
    
# define the derivative of f(x)
def f_prime(x):
    return -((A / deltatime) * x) - ((W*(S)**0.5)/n) * (((W*x)/(W + 2*x))**0.66 + (0.66 * x) * (((1/(x)) + (2/W))**0.33) * (W**2)/((W+(2*x))**2))

# define the initial guess
x0 = 101
# define the tolerance
tol = 1e-6
# define the maximum number of iterations
max_iter = 100
# initialize the iteration counter
n = 0
# compute the initial error
error = abs(f(x0))

# iterate until the error is less than the tolerance
# or the maximum number of iterations is reached
while error > tol and n < max_iter:
    # update the approximation
    x1 = x0 - f(x0) / f_prime(x0)

    # update the error
    error = abs(f(x1))

    # update the iteration counter
    n += 1

    # update the initial guess
    x0 = x1

# print the final approximation
print(x1)

错误定位 #

错误触发行：

def f(x):
    return qin_secs-(((x-b)/deltatime)*A)-((W*x)/n)*((W*x)/((W+2)*x))**0.66*(S)**0.5

尝试的无效修复：

def f(x):
    if x == 0:
        return 0
    return qin_secs-(((x-b)/deltatime)*A)-((W*x)/n)*((W*x)/((W+2)*x))**0.66*(S)**0.5

错误根源与修复方案 #

1. 核心问题：变量名冲突 #

全局变量中定义了曼宁系数n = 0.04，但后续迭代时又将迭代计数器命名为n并赋值为0，导致f(x)函数中(W*x)/n变为除以0，触发除零错误——和x的值无关。

2. 修复步骤 #

• 修改迭代计数器变量名：将迭代计数器n改为iter_count，避免与全局变量n重名
• 简化f(x)表达式：(W*x)/((W+2)*x)在x≠0时可约分为W/(W+2)，既简化计算又避免x相关的潜在除零问题
• 修正x=0时的返回值：原修复中返回0不符合函数逻辑，应代入x=0计算正确结果

修改后的完整代码 #

# Newton-Raphson iterative scheme
qin_secs = 100 
W = 3
L = 7000
A = W * L
n = 0.04  # 曼宁系数，保留原变量名
S = 0.007
b = 100
deltatime = 1

# 定义目标函数f(x)
def f(x):
    if x == 0:
        # 代入x=0计算正确值，而非返回0
        return qin_secs + (b * A) / deltatime
    # 简化表达式：约去x，避免不必要的计算和潜在错误
    flow_term = (W / n) * (W / (W + 2))**0.66 * (S)**0.5 * x
    return qin_secs - ((x - b) * A) / deltatime - flow_term
    
# 定义导数函数f_prime(x)
def f_prime(x):
    if x == 0:
        # 处理x=0时的导数，避免除零
        return -A / deltatime - (W*(S)**0.5)/n * (0 + 0)
    # 简化导数中的表达式
    ratio = W * x / (W + 2 * x)
    term1 = ratio**0.66
    term2 = 0.66 * x * ((1/x + 2/W)**0.33) * (W**2) / ((W + 2*x)**2)
    return -A / deltatime - ((W*(S)**0.5)/n) * (term1 + term2)

# 初始猜测值
x0 = 101
# 容差
tol = 1e-6
# 最大迭代次数
max_iter = 100
# 初始化迭代计数器，避免变量名冲突
iter_count = 0
# 初始误差
error = abs(f(x0))

# 迭代循环
while error > tol and iter_count < max_iter:
    x1 = x0 - f(x0) / f_prime(x0)
    error = abs(f(x1))
    iter_count += 1
    x0 = x1

# 输出结果
print(f"迭代结果：{x1}")

内容的提问来源于stack exchange，提问作者dinn_