能否用Python求解该方程组？请求代码正确性验证

我尝试用Python编写方程组求解代码，但几乎没有相关经验。已写出代码，但不确定其正确性，希望得到帮助。待求解的方程组对应参考图片。

我的代码如下：

import numpy as np
from scipy.optimize import fsolve

# Initial data
Ptmax_init = 100
Q_init = 100
e_init = 0.01

r0 = 0.85
S = 0.18
b = 0.045
c0 = 9.8 * 10**3
B = 0.53
th = 0.05
Cr = 6.33 * 10**5
mu = 0.8
Qc = 1.177 * 10**4
h_kr = 0.6
L = 2.5
ec = 0.01
k = 2.81
tang = 0.65

#Equations
def equations(vars):
    Ptmax, Q, e = vars
    eq1 = Ptmax = Q * ((r0 / (r0 - e)) * ((S - b) / S) * ((2 * c0 * (B + 2 * th) * r0) / (Cr * e) + tang - ((2 * c0 * (B + 2 * th)) / Cr)) + (mu * (b / S)))
    eq2 = Q = (Qc + (Ptmax * h_kr / L))
    eq3 = e = (ec + (Ptmax * h_kr) / (L * 1.44 * np.sqrt((Cr**2 * Q) / r0 * (1 + np.sqrt(k + 1))**2)**(1/3)))
    return [eq1, eq2, eq3]

# Solving
initial_guess = [Ptmax_init, Q_init, e_init]
solution = fsolve(equations, initial_guess, maxfev=150)
Ptmax_sol, Q_sol, e_sol = solution

print(f"Solve: Ptmax = {Ptmax_sol}, Q = {Q_sol}, e = {e_sol}")

代码问题分析与修正建议 #

1. 方程组定义错误
   fsolve要求传入的方程组是等于0的等式，但你的代码里用了赋值操作（eq1 = Ptmax = ...），这会导致返回的是赋值后的变量值，而非方程的残差。正确写法是将右边的表达式移到左边，让方程等于0：

   • eq1 = Ptmax - Q * ((r0/(r0-e)) * ((S-b)/S) * ((2*c0*(B+2*th)*r0)/(Cr*e) + tang - (2*c0*(B+2*th))/Cr) + mu*(b/S))
   • eq2 = Q - (Qc + Ptmax*h_kr/L)
   • eq3 = e - (ec + (Ptmax*h_kr)/(L * 1.44 * ((Cr**2 * Q / r0) * (1 + np.sqrt(k+1))**2)**(1/6)))

2. 运算优先级优化
   eq3中的np.sqrt(...)**(1/3)可以简化为(...)**(1/6)（因为平方根是1/2次方，再开三次方就是1/2 * 1/3 = 1/6次方），既避免嵌套根号的优先级问题，也提升计算效率。

3. 初始猜测值合理性调整
   原初始猜测Q_init=100远小于Qc=1.177e4，会导致fsolve难以收敛。建议将初始猜测调整为更接近真实解的数值：

   • initial_guess = [100, 1.2e4, 0.01]（Q的初始值接近Qc，Ptmax根据eq2估算）

修正后的代码 #

import numpy as np
from scipy.optimize import fsolve

# 初始参数
Ptmax_init = 100
Q_init = 1.2 * 10**4
e_init = 0.01

r0 = 0.85
S = 0.18
b = 0.045
c0 = 9.8 * 10**3
B = 0.53
th = 0.05
Cr = 6.33 * 10**5
mu = 0.8
Qc = 1.177 * 10**4
h_kr = 0.6
L = 2.5
ec = 0.01
k = 2.81
tang = 0.65

# 定义方程组（残差形式，即方程=0）
def equations(vars):
    Ptmax, Q, e = vars
    # 方程1：Ptmax - 右边表达式 = 0
    eq1 = Ptmax - Q * (
        (r0 / (r0 - e)) * ((S - b) / S) * 
        ((2 * c0 * (B + 2 * th) * r0) / (Cr * e) + tang - (2 * c0 * (B + 2 * th)) / Cr) + 
        mu * (b / S)
    )
    # 方程2：Q - 右边表达式 = 0
    eq2 = Q - (Qc + Ptmax * h_kr / L)
    # 方程3：e - 右边表达式 = 0，简化根号运算
    term = ((Cr**2 * Q / r0) * (1 + np.sqrt(k + 1))**2)**(1/6)
    eq3 = e - (ec + (Ptmax * h_kr) / (L * 1.44 * term))
    return [eq1, eq2, eq3]

# 求解方程组
initial_guess = [Ptmax_init, Q_init, e_init]
solution = fsolve(equations, initial_guess)
Ptmax_sol, Q_sol, e_sol = solution

print(f"求解结果：Ptmax = {Ptmax_sol:.2f}, Q = {Q_sol:.2f}, e = {e_sol:.6f}")

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