如何使用fsolve求解并绘制4-8kg食材的烹饪时间（小时）曲线？

数学模型梳理 #

先明确核心方程（原无穷级数取k=22项）：

57/80 = (2/(0.464×π)) × Σₖ=1²² [ (-1)^(k+1)/k × sin(kπ×0.464) × exp( -α×(kπ/a)²×t ) ]

参数定义：

• weight：4~8kg的数组
• a = (weight/(1360×π))^(1/3)（单位：m，默认weight单位为kg，1360为密度kg/m³）
• α = 0.7×10^-7（热扩散系数，单位：m²/s）
• 目标：对每个weight求解对应烹饪时间t（秒），转换为小时后绘制曲线。

Python 实现代码 #

import numpy as np
from scipy.optimize import root_scalar
import matplotlib.pyplot as plt

# 固定参数
ALPHA = 0.7e-7
TARGET = 57/80
R_OVER_A = 0.464
COEFF = 2 / (R_OVER_A * np.pi)
K_MAX = 22

# 生成weight数组（4到8kg，取50个点保证曲线平滑）
weights = np.linspace(4, 8, 50)
cook_time_hours = np.zeros_like(weights)

# 定义每个weight对应的方程：f(t) = 右边计算值 - 目标值，求f(t)=0的根
def equation(t, a):
    k = np.arange(1, K_MAX+1)
    # 向量化计算求和项，避免循环维度错误
    sum_term = np.sum(
        ((-1)**(k+1)/k) * np.sin(k * np.pi * R_OVER_A) * 
        np.exp(-ALPHA * (k * np.pi / a)**2 * t)
    )
    return COEFF * sum_term - TARGET

# 遍历每个weight求解对应烹饪时间
for i, w in enumerate(weights):
    a = (w / (1360 * np.pi)) ** (1/3)
    # 用brentq方法求解，设置合理区间避免发散（1小时到24小时）
    sol = root_scalar(equation, args=(a,), method='brentq', bracket=[3600, 86400])
    if sol.converged:
        cook_time_hours[i] = sol.root / 3600  # 转换为小时
    else:
        cook_time_hours[i] = np.nan  # 标记未收敛的异常点

# 绘制曲线
plt.figure(figsize=(10,6))
plt.plot(weights, cook_time_hours, 'b-', linewidth=2)
plt.xlabel('Weight (kg)')
plt.ylabel('Cooking Time (hours)')
plt.title('Cooking Time vs Meat Weight')
plt.grid(True)
plt.show()

关键问题修复说明 #

1. 数组维度匹配：
   • 对每个weight单独计算a，再传入单变量方程求解，避免广播维度不匹配问题。
   • 用np.arange生成k数组，向量化计算求和项，替代循环计算，既高效又避免维度错误。
2. 求解稳定性：
   • 使用brentq区间求解法，传入合理的时间区间（1-24小时），比单一初始猜测更稳定，避免求解发散。
3. 结果处理：
   • 统一将求解得到的秒数转换为小时，符合绘图的单位需求；标记未收敛的异常点，便于排查数据问题。

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