多元Newton法求解四元非线性方程组Python程序无输出问题求助

代码核心错误点 #

• 雅可比矩阵计算逻辑完全错误：autograd.jacobian的输入是待求导的函数，不是函数计算得到的标量结果，你当前代码中将f1~f4的计算值直接传入jacobian，根本无法得到正确的导数，且强行将4个导数结果reshape为4x1矩阵，既不符合4维方程组雅可比矩阵应为4x4的要求，也会导致后续求逆、求解运算直接报错或卡死。
• 牛顿迭代公式参数顺序完全颠倒：np.linalg.solve(A, b)的作用是求解线性方程组Ax = b，牛顿法标准更新逻辑为x_new = x_old - J^{-1} * F(x_old)，等价于求解J * delta = -F(x_old)，你当前代码将函数值矩阵和雅可比矩阵的传入顺序完全搞反，逻辑不成立。
• 收敛判断逻辑顺序错误：你在计算两次迭代的欧式距离前，已经将x0赋值为新的迭代值newxvalues，导致newxvalues - x0永远为0，循环判断逻辑完全失效。
• 无任何结果输出逻辑：迭代完成后没有打印任何结果，自然看不到运行输出。

修复后可运行代码 #

from autograd import jacobian
import autograd.numpy as np    
import numpy.linalg 
import sys 

sys.setrecursionlimit(10000)

def functionmatrix(x0):
    f1 = pow((x0[0] - 1560), 2) + pow((x0[1] - 6540), 2) + pow((x0[2] - 20140), 2) - (2.9*10**3*pow((0.07074 - x0[3]), 2))
    f2 = pow((x0[0] - 18760), 2) + pow((x0[1] - 2750), 2) + pow((x0[2] - 18610), 2) - (2.9*10**3*pow((0.07220 - x0[3]), 2))
    f3 = pow((x0[0] - 17610), 2) + pow((x0[1] - 14630), 2) + pow((x0[2] - 13480), 2) - (2.9*10**3*pow((0.07690 - x0[3]), 2))
    f4 = pow((x0[0] - 19170), 2) + pow((x0[1] - 610), 2) + pow((x0[2] - 18390), 2) - (2.9*10**3*pow((0.07242 - x0[3]), 2))
    return np.array([f1, f2, f3, f4]).reshape(4,1)

# 直接通过autograd生成雅可比计算函数，无需手动重复实现
jac = jacobian(functionmatrix)

def determinevalues():
    euclideandistance = 100
    iteration = 0
    tolerance = 1e-8
    maxiterations = 1000
    x0 = np.array([0,0,6370,0], dtype = float).reshape(4,1)
    while(np.any(abs(euclideandistance) > tolerance) and iteration < maxiterations):
        F = functionmatrix(x0)
        J = jac(x0).reshape(4,4) # 生成4x4雅可比矩阵
        delta = np.linalg.solve(J, -F) # 求解迭代增量
        newxvalues = x0 + delta
        euclideandistance = np.linalg.norm(newxvalues - x0)
        iteration += 1
        x0 = newxvalues
        # 可选：打印迭代过程
        print(f"迭代次数：{iteration}，残差范数：{euclideandistance:.6f}")
    # 输出最终结果
    print("\n求解结果：")
    print(f"x = {x0[0][0]:.4f}")
    print(f"y = {x0[1][0]:.4f}")
    print(f"z = {x0[2][0]:.4f}")
    print(f"d = {x0[3][0]:.4f}")
    return x0

def main():
    determinevalues()
    
if __name__ == "__main__":
    main()

额外优化说明 #

• 直接使用autograd提供的雅可比生成接口，无需重复编写函数逻辑，减少出错概率
• 不用显式计算雅可比矩阵的逆，用np.linalg.solve直接求解增量，数值稳定性更高
• 增加了迭代过程打印，方便观察收敛情况
• 调整了迭代逻辑顺序，收敛判断完全符合要求

内容的提问来源于stack exchange，提问作者Weightlifting Without Limits