Python实现Newton法迭代异常问题排查

原代码 #

p= lambda x: 3 * x**3 - 5.03 * x**2 - 1.95 * x + 0.02

dp = lambda x: 9 * x**2 - 10.06 * x - 1.95
x0 = -1
tol = 10**-4
max_iter = 100

def newton_method(f, df, x0, tol, max_iter=100):
    for i in range(max_iter):
        fx = f(x0)
        dfx = df(x0)
        if df==0 and dfx!=0 :
            raise Exception('division by 0')
        
        elif fx==0:
            return i, x0
        
        x = x0 - fx/dfx
        
        if abs(x-x0) <=tol:
              return i, x0
        x0 = x
    raise ValueError("The newton method doesn't converge after {} iterations".format(max_iter))

c,r = newton_method(p, dp , x0, tol, max_iter)
print("Approssimazione:", r,c)

问题分析与错误点 #

1. 导数为零的判断逻辑完全错误 #

原代码中if df==0 and dfx!=0的逻辑完全颠倒：df是传入的导数函数，永远不可能等于0；正确逻辑应该是判断导数值dfx是否接近0（浮点运算中不能直接用==0），避免除以极小值导致迭代发散。

2. 函数值为零的判断不严谨 #

直接用fx==0在浮点运算中几乎无法触发（除非刚好是精确零点），应该通过判断函数值的绝对值是否小于极小阈值（比如1e-12）来识别接近零点的情况。

3. 收敛判断返回错误值 #

原代码在满足abs(x-x0) <= tol时返回旧值x0，但x才是迭代后的新近似解，导致返回结果并非最终迭代的精度值；同时迭代次数计数有误，第一次迭代完成后应计数为1而非0。

4. 初始点5附近收敛到2的原因 #

首先明确：x=2是函数的精确零点（手动计算验证：3*2³ -5.03*2² -1.95*2 +0.02 = 24 -20.12 -3.9 +0.02 = 0）。初始点5附近收敛到2是牛顿法的特性——迭代收敛到哪个根取决于初始点所在的吸引域，若要收敛到5.3附近的根，需将初始点调整至更靠近5.3的位置（如x0=5.2），或使用带阻尼的修正牛顿法。

修正后的代码 #

p = lambda x: 3 * x**3 - 5.03 * x**2 - 1.95 * x + 0.02
dp = lambda x: 9 * x**2 - 10.06 * x - 1.95

def newton_method(f, df, x0, tol, max_iter=100):
    for i in range(max_iter):
        fx = f(x0)
        dfx = df(x0)
        
        # 判断导数值是否接近0，避免除以极小值
        if abs(dfx) < 1e-12:
            raise Exception('导数接近0，无法继续迭代')
        
        # 判断当前点是否足够接近零点
        if abs(fx) < tol:
            return i, x0
        
        # 牛顿迭代计算下一个点
        x = x0 - fx / dfx
        
        # 判断迭代步长是否满足精度要求，返回新的近似解
        if abs(x - x0) <= tol:
            return i+1, x
        x0 = x
    
    raise ValueError(f"牛顿法在{max_iter}次迭代后未收敛")

# 测试初始点x0=5.2（靠近5.3的根）
c, r = newton_method(p, dp, 5.2, 1e-4, 100)
print(f"近似解：{r}，迭代次数：{c}")

# 测试初始点x0=2
c2, r2 = newton_method(p, dp, 2, 1e-4, 100)
print(f"近似解：{r2}，迭代次数：{c2}")

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