The problem I am facing consists of plotting the graph I attach, only the discontinuous lines at the moment.
.
For that purpose I need to take into account these equations:
.
Before starting explaining the process, there is one constant, gamma=1.4, which I am going to be using.
To start plotting, I take equation 1. I take values of M(L) from 0.1 to 1 and of "x" (x=2fL/De) from 10e-1 to 10e3 (in a continuous way), as you can see in the graph. From here all I do is get the value of M(0) for each of the values of M(L), so that I can have 10 different curves later.
Once I have all M(0) values, I go to the second equation and by introducing the values of M(0) and M(L) I obtain the corresponding values of the pressure ratio (P0(0)/P(L))--> y axis.
For a reason I cannot explain, when asking Chatgpt for the python code, I get the following one copied below, with its corresponding graph, and the curves do not match with the ones in the image 1. I have extracted directly the latex code with Mathpix from the equations involved for chatgpt to convert it into python code, and that way i assure that it is writing the equations right.
Can someone help me get to the correct resolution of my problem?
The code:
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()
I then tried going step by step, calculating first the relationship between M(0) and x. I cleared the value of M(0) as in image 2 [image 2][1] [1]: https://i.stack.imgur.com/D0pwT.png
I got the code performing the numerical calculation with Jacobi method, as follows:
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
