压力比与摩擦力和其他几何因素

问：

我面临的问题包括绘制我所附的图形，目前只有不连续的线。

.

为此，我需要考虑以下方程式：

.

在开始解释这个过程之前，有一个常数，gamma=1.4，我将使用它。

为了开始绘图，我取等式 1。我从 0.1 到 1 取 M（L） 的值，从 10e-1 到 10e3 取“x”（x=2fL/De） 的值（以连续的方式），如图所示。从这里开始，我所要做的就是获得 M（L） 的每个值的 M（0） 值，以便以后可以有 10 条不同的曲线。

一旦我得到了所有 M（0） 值，我就进入第二个方程，通过引入 M（0） 和 M（L） 的值，我得到了压力比 （P0（0）/P（L））--> y 轴的相应值。

由于我无法解释的原因，当向 Chatgpt 询问 python 代码时，我得到了下面复制的以下代码及其相应的图表，并且曲线与图 1 中的曲线不匹配。我使用 Mathpix 直接从 chatgpt 所涉及的方程式中提取了 latex 代码，以将其转换为 python 代码，这样我就可以确保它正确地编写方程式。

有人可以帮助我正确解决我的问题吗？

代码：

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()

然后我试着一步一步地走，首先计算 M（0） 和 x 之间的关系。我清除了M（0）的值，如图2所示[图2][1] [1]：https://i.stack.imgur.com/D0pwT.png

我得到了使用 Jacobi 方法执行数值计算的代码，如下所示：

import matplotlib.pyplot as plt

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

    M0_prev = 1.0  # Valor inicial para M(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()

I am getting some errors regarding "overflow encountered in scalar multiply" and "invalid value encountered in scalar divide" plus some of the results I get for M0 are NaN and there should be 1000 values for M0 since there are 1000 values of x, which is not my case. Im pretty sure that my equation is right, but maybe I am not using the method properly