求助：压力比与摩擦及几何参数曲线绘制的代码修正

• 需绘制压力比**P₀(0)/P(L)随2fL/De（记为x）**变化的虚线曲线，已知常数γ=1.4
• 计算逻辑：
  1. 取M(L)（0.1至1）和x（1e-1至1e3）的连续值，求解M(0)
  2. 将M(0)代入对应方程得到压力比值
• 遇到的问题：
  • ChatGPT生成的代码绘制曲线与目标不符
  • 自行编写的雅可比求解代码出现overflow encountered in scalar multiply、invalid value encountered in scalar divide错误，部分M(0)结果为NaN，且无法生成与x数量匹配的1000个M(0)值，确认方程正确但求解方法使用不当

import numpy as np
import matplotlib.pyplot as plt

# Function to solve for M(0) given M(L)
def solve_M0(M_L):
    y = np.log((x**2 / M_L**2) * (1 + (gamma - 1) * M_L**2 / 2) /
               (1 + (gamma - 1) * x**2 / 2)) + 2 / (gamma + 1) * (1 / x**2 - 1 / M_L**2)
    return np.exp(y)
    
# Function to solve for p0(0)/p(L) given M(0)
def solve_p_ratio(M_0):
    y = (M_0**2 / (1 + (gamma - 1) / 2 * M_0**2)**((gamma + 1) / (gamma - 1))) * x**2 / (1 + (gamma - 1) / 2 * x**2)
    return np.sqrt(y)

# Plotting
M_L_values = np.arange(0.1, 1.1, 0.1)  # Values of M_L
x = np.linspace(1e-1, 1e3, 1000)  # Values for x-axis (2*f*L/De)
# Constants
gamma = 1.4

plt.figure(figsize=(10, 6))
for M_L in M_L_values:
    M_0 = solve_M0(M_L)
    p_ratio = solve_p_ratio(M_0)
    plt.plot(x, p_ratio, label=f'M_L = {M_L}')

plt.xlabel('2*f*L/De')
plt.ylabel('p0(0)/p(L)')
plt.xscale('log')
plt.legend()
plt.grid(True)
plt.title('Plot of p0(0)/p(L) vs 2*f*L/De for different M_L values')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def jacobi_method(ML, x):
    gamma = 1.4
    tol = 1e-6
    max_iter = 100

    M0_prev = 1.0  # 初始值设为1.0
    M0 = M0_prev

    for _ in range(max_iter):
        M0 = np.sqrt(ML **2 * (1 + (gamma - 1) / 2 * M0 **2 * 
        np.exp(gamma / (gamma + 1) * x) * np.exp(-2 * (ML**2 - M0**2) / 
        (M0**2 * ML**2 * (gamma + 1))))) / (1 + (gamma - 1) / 2 * ML**2)
        
        if np.abs(M0 - M0_prev) < tol:
            break

        M0_prev = M0

    return M0

ML_values = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
x_values = np.logspace(-1, 2, num=1000)

M0_values = np.zeros((len(ML_values), len(x_values)))

for i, ML in enumerate(ML_values):
    for j, x in enumerate(x_values):
        M0_values[i, j] = jacobi_method(ML, x)

plt.figure(figsize=(10, 6))
plt.xlabel('x')
plt.ylabel('M(0)')

for i, ML in enumerate(ML_values):
    plt.plot(x_values, M0_values[i, :], label=f'M(L) = {ML}')

plt.legend()
plt.grid(True)
plt.show()

报错信息 #

• overflow encountered in scalar multiply
• invalid value encountered in scalar divide
• 部分M(0)结果为NaN，无法生成与x数量匹配的1000个有效结果

问题分析 #

1. ChatGPT代码问题：solve_M0函数错误地将x作为全局变量使用，且推导的方程不符合管道流动中M(0)与M(L)的物理关系（正确模型应为Fanno流摩擦公式）。
2. 雅可比方法问题：
   • 迭代公式推导有误，指数项易引发数值溢出
   • 初始值选择不当（亚声速流中M(0)应小于M(L)，设为1.0会导致迭代发散）
   • 未处理除零、溢出等异常情况

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt

gamma = 1.4

# Fanno流标准摩擦函数：给定马赫数M，计算对应的Fanno参数值
def fanno_function(M):
    term1 = (gamma + 1)/(2*gamma) * np.log( ((gamma + 1)*M**2) / (2 + (gamma - 1)*M**2) )
    term2 = (1 - M**2)/(gamma*M**2)
    return term1 + term2

# 二分法求解M(0)：已知M(L)和x=2fL/De，满足fanno(M0) - fanno(ML) = x
def solve_M0(ML, x):
    # 亚声速流中M0 < ML < 1，设置合理搜索区间
    low = 1e-6
    high = ML
    tol = 1e-8
    max_iter = 100
    
    for _ in range(max_iter):
        mid = (low + high)/2
        f_mid = fanno_function(mid) - fanno_function(ML) - x
        f_low = fanno_function(low) - fanno_function(ML) - x
        
        if np.sign(f_mid) == np.sign(f_low):
            low = mid
        else:
            high = mid
            
        if np.abs(f_mid) < tol:
            break
    
    return (low + high)/2

# 计算压力比P0(0)/P(L)
def compute_pressure_ratio(M0, ML):
    # 总压比P0(0)/P0(L)
    p0_ratio = ( ( (gamma+1)*ML**2 / (2 + (gamma-1)*ML**2) ) / ( (gamma+1)*M0**2 / (2 + (gamma-1)*M0**2) ) ) ** (gamma/(gamma-1))
    # 静总压比P0(L)/P(L)
    p0L_pL = (1 + (gamma-1)/2 * ML**2)**(gamma/(gamma-1))
    # 最终压力比P0(0)/P(L)
    return p0_ratio * p0L_pL

# 生成计算参数
ML_values = np.arange(0.1, 1.0, 0.1)
x_values = np.logspace(-1, 3, 1000)

# 绘制虚线曲线
plt.figure(figsize=(10, 6))
for ML in ML_values:
    p_ratios = []
    for x in x_values:
        M0 = solve_M0(ML, x)
        p_ratio = compute_pressure_ratio(M0, ML)
        p_ratios.append(p_ratio)
    plt.plot(x_values, p_ratios, linestyle='--', label=f'M(L) = {ML:.1f}')

plt.xlabel('2fL/De')
plt.ylabel('P₀(0)/P(L)')
plt.xscale('log')
plt.legend()
plt.grid(True, which='both', linestyle='--', alpha=0.7)
plt.title('压力比P₀(0)/P(L)随2fL/De变化曲线')
plt.show()

修正说明 #

1. 采用Fanno流标准摩擦公式构建物理模型，确保计算逻辑符合流体力学规律
2. 使用二分法求解M(0)，替代不稳定的雅可比方法，避免迭代发散和数值溢出
3. 明确压力比的推导逻辑，结合总压比与静总压比计算最终结果
4. 设置合理的搜索区间与收敛阈值，避免NaN和无效结果
5. 绘制虚线曲线，满足需求格式要求

内容的提问来源于stack exchange，提问作者Néstor Indu