如何在Python中用牛顿法求解方程组？兼容单变量场景

你的核心需求是让牛顿法同时支持单变量方程和固定部分变量的多元方程组求解，不需要把方程组转为单个方程，只需要重构原函数，增加变量固定逻辑并兼容向量/矩阵运算。

原代码的核心问题 #

1. 依赖eval传入函数字符串，既不安全也不灵活，无法直接处理返回向量的方程组
2. 仅支持标量除法，无法处理多元方程组需要的雅可比矩阵求逆/伪逆运算
3. 没有变量固定的逻辑，无法实现“固定三个变量、求解单个变量”的需求

修改后的通用牛顿法函数 #

使用numpy处理向量/矩阵运算，去掉eval直接传入函数对象，增加固定变量的参数：

import numpy as np

def fnewton(f, df, x0, n_iter, fixed_vars=None):
    """
    通用牛顿法，支持单变量和固定部分变量的多元方程组求解
    参数：
        f: 目标函数，输入向量x，返回标量（单变量）或向量（多元残差）
        df: 导数/雅可比矩阵函数，输入向量x，返回标量（导数）或矩阵（雅可比）
        x0: 初始值向量
        n_iter: 迭代次数
        fixed_vars: 字典，键为固定变量的索引，值为固定值（如{1:0.5,2:100}表示固定x[1],x[2]）
    """
    x = np.array(x0, dtype=np.float64)
    fixed_idxs = []
    fixed_vals = []
    
    if fixed_vars:
        fixed_idxs = list(fixed_vars.keys())
        fixed_vals = list(fixed_vars.values())
    
    for _ in range(n_iter):
        # 每次迭代前恢复固定变量的值
        for idx, val in zip(fixed_idxs, fixed_vals):
            x[idx] = val
        
        res = np.array(f(x))
        jac = np.array(df(x))
        
        # 自动区分单变量和多元场景
        if res.ndim == 0 and jac.ndim == 0:
            # 单变量标量逻辑
            delta = res / jac
            x = x - delta
        else:
            # 多元场景用最小二乘伪逆，避免矩阵不可逆
            delta = np.linalg.lstsq(jac, res, rcond=None)[0]
            x = x - delta
    
    print('根的近似值：')
    print(x)
    print(f'迭代次数：{n_iter}')
    return x

针对你的方程组的使用示例 #

假设你要固定y、nl、nv，求解T（x[0]），需要先实现方程组的雅可比矩阵函数：

1. 实现雅可比矩阵函数 #

def dc_dx(x):
    T = x[0]
    y = x[1]
    nl = x[2]
    nv = x[3]
    
    # 实现Antoine函数对T的导数（根据Antoine公式推导）
    # Antoine公式一般为 log10(P) = A - B/(T+C)，导数dP/dT = P * ln(10) * B/(T+C)^2
    def dAntoine_dT(A, B, C, T):
        P = 10 ** (A - B/(T + C))
        return P * np.log(10) * B / (T + C)**2
    
    dAntoineN_dT = dAntoine_dT(An, Bn, Cn, T)
    dAntoineH_dT = dAntoine_dT(Ah, Bh, Ch, T)
    
    # 构造雅可比矩阵：每行对应一个方程对四个变量的偏导数
    jac = [
        [0.63 * dAntoineN_dT, -760, 0, 0],    # RLN对T、y、nl、nv的偏导
        [(1-0.63)*dAntoineH_dT, 760, 0, 0],   # RLH对T、y、nl、nv的偏导
        [0, nv, 0.63, y],                     # N对T、y、nl、nv的偏导
        [0, -nv, (1-0.63), (1-y)]             # H对T、y、nl、nv的偏导
    ]
    return jac

2. 调用牛顿法求解 #

# 初始值：[T初始值, y固定值, nl固定值, nv固定值]
initial_x = [300, 0.5, 100, 200]
# 指定固定变量：索引1(y)=0.5，索引2(nl)=100，索引3(nv)=200
fixed_vars = {1: 0.5, 2: 100, 3: 200}

# 调用函数求解T
result = fnewton(c, dc_dx, initial_x, 10, fixed_vars=fixed_vars)

关键说明 #

• 如果只需要让方程组中的某一个方程满足残差为0，可以修改c函数，只返回该方程的结果，此时会自动退化为单变量求解逻辑
• 用np.linalg.lstsq代替直接求逆，能避免雅可比矩阵不可逆的情况，提升鲁棒性
• 去掉eval后，代码安全性和可维护性大幅提升，也能直接处理任意复杂度的函数

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