如何使用scipy求解给定A_e/A*列表对应M_e的非线性方程？

原有代码错误点 #

• 目标函数返回值错误：result.any()返回布尔值（数组存在非零值即返回True，对应数值1），scipy.optimize.newton要求返回残差的实数值，布尔值隐式转0/1后完全无法支撑收敛判断。
• 批量输入处理逻辑错误：newton为单根求解器，每次仅能求解单个方程的根，将整个A_e/A*数组写入目标函数相当于同时求解数十个独立方程，必然无法收敛。
• 公式指数计算错误：喷管截面积比公式的指数应为(gamma+1)/(2*(gamma-1))，原代码误写为(gamma+1)/(2*gamma-1)，公式本身错误导致残差计算完全偏离正确值。

修正后代码 #

import numpy as np
from scipy.optimize import newton

# 固定参数
gamma = 1.2
# 生成A_e/A*输入列表
A_ratio_list = np.arange(1, 1.25, step=0.004)
M_e_result = []

# 定义单值残差函数，x为马赫数，A_ratio为对应输入的截面积比
def residual(x, A_ratio):
    inner_term = (2 / (gamma + 1)) * (1 + (gamma - 1) / 2 * x ** 2)
    calc_A_ratio = (1 / x) * (inner_term) ** ((gamma + 1) / (2 * (gamma - 1)))
    return calc_A_ratio - A_ratio

# 遍历每个输入值单独求解
for A_r in A_ratio_list:
    # 初始 guess 设为1.1求解超声速解，若需要亚声速解可改为0.5
    me = newton(residual, 1.1, args=(A_r,))
    M_e_result.append(me)

# 转为numpy数组方便后续计算
M_e_result = np.array(M_e_result)

结果验证 #

你可以通过如下代码打印对应关系验证正确性：

for a_r, me in zip(A_ratio_list, M_e_result):
    print(f"A_e/A* = {a_r:.3f}，对应M_e = {me:.4f}")

如果需要更高的求解稳定性，可替换newton为scipy.optimize.root_scalar，指定求解边界避免发散。

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