Python实现双方程牛顿法：Julia求根代码转写求助

逻辑说明 #

这段代码不存在双方程求解需求，属于典型的单变量单方程求根+公式代入计算场景：

• 内部定义的β_func是嵌套的单变量函数，仅需要在指定区间内找到让函数值为0的根β即可
• Julia代码中find_zero(β_func, (0, 0.1), Bisection())调用的是区间二分法求根，属于最基础的单变量求根方法，不需要参考多方程联立求解的案例
• 求出β后代入后续算术公式就能得到最终结果

前置准备 #

首先安装数值计算依赖，执行命令：
pip install numpy scipy

完整实现代码 #

代码逻辑和原Julia代码逐行对应，仅对x=0处的除零问题做了数值兼容处理：

import numpy as np
from scipy.optimize import root_scalar

# 以下常量替换为原Julia脚本中对应的实际参数值即可
P_1 = None
T_AP = None
T_P1 = None
T_APL = None
T_PL = None
T_APH = None
T_PH = None
alpha = None
I_sat = None

def p(temp, irr):
    temp_k = temp + 273.15
    # 计算温度响应的最大光合速率
    p_max = P_1 * np.exp(T_AP / T_P1 - T_AP / temp_k) / (
        1 
        + np.exp(T_APL / temp_k - T_APL / T_PL) 
        + np.exp(T_APH / T_PH - T_APH / temp_k)
    )

    # 定义求根目标函数
    def beta_func(x):
        calc_val = (alpha * I_sat / np.log(1 + alpha / x)) * (alpha / (alpha + x)) * ((x / (alpha + x)) ** (x / alpha))
        return p_max - calc_val

    # 二分法求β，左端点取极接近0的正数避免x=0除零报错
    root_res = root_scalar(beta_func, bracket=[1e-10, 0.1], method="bisect")
    beta = root_res.root

    # 计算最终光合值
    p_s = alpha * I_sat / np.log(1 + alpha / beta)
    return p_s * (1 - np.exp(-alpha * irr / p_s)) * np.exp(-beta * irr / p_s)

常见问题说明 #

• 原Julia代码给出的求根区间是(0, 0.1)，但x=0时alpha/x会触发除零错误，取1e-10作为左端点不会影响求根精度，和原代码运行结果差异可以忽略
• 如果需要批量计算多组温度、辐照对应的结果，直接传入numpy数组格式的temp、irr即可，numpy广播机制会自动完成批量计算，不需要手写循环
• 如果运行时提示二分法区间端点函数值同号，先检查你填入的常量参数是否和原Julia代码一致，正常情况下两个端点的函数值符号相反，满足二分法的使用条件

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