如何用Python绘制雷诺数与剪切应力的双对数散点图

Hey there! I see where you're stuck—let's break this down and get your 100,000-point log-log scatter plot working properly.

First, let's recap the gaps in your current code:

• You're only generating 10 shear stress values per Reynolds Number (Re) instead of 10,000
• There's no logic to collect all Re and shear stress pairs into a format that can be plotted
• You haven't added code to render the log-log scale plot itself

Here's a revised, optimized solution using NumPy and Matplotlib to handle the large dataset efficiently:

Step 1: Generate 100 Reynolds Numbers (0.1 to 2017) #

Instead of a for-loop with 0.1*(1.1**i), we use np.geomspace to create exactly 100 evenly spaced values on a logarithmic scale between 0.1 and 2017. This ensures we hit your desired range perfectly.

Step 2: Generate 10,000 Shear Stress Values per Re #

We leverage vectorized NumPy operations instead of nested for-loops—this is way faster for large datasets. For each Re, we'll generate 10,000 shear stress values in one go.

Step 3: Collect Data for Plotting #

We repeat each Re 10,000 times to match the number of shear stress values, so every data point has a corresponding x (Re) and y (shear stress) value.

Full Code #

import numpy as np
import matplotlib.pyplot as plt
import math

# Replace these placeholders with your actual constant values!
mean = 0  # Your mean for normal distribution
sd = 1    # Your standard deviation
fi = 30   # Your fi value (converted to radians below if in degrees)
d50 = 0.1 # Your d50 parameter
D = 1.0   # Your D parameter

# Step 1: Generate 100 logarithmically spaced Reynolds Numbers
Re_values = np.geomspace(0.1, 2017, num=100)

# Initialize lists to store all data points
all_Re = []
all_shear_stress = []

# Step 2: Generate 10,000 shear stress values for each Re
for Re in Re_values:
    # Precompute reusable terms for this Re
    exp_term_08 = np.exp(-0.08 * Re)
    exp_term_10 = np.exp(-0.1 * Re)
    
    B = exp_term_08 * (2.5 * np.log(Re) + 5.25) + 8.5 * (1 - exp_term_08)
    C = 0.8 + 0.9 * ((exp_term_08 / (Re**2)) + ((1 - exp_term_08) / (B**2))) ** (-0.5)
    F = 0.31 * Re * exp_term_10 + 1.8 * np.exp(-0.88 * d50 / D) * (1 - exp_term_10)
    A = F / C
    
    # Generate 10,000 normal samples (vectorized, no loop needed!)
    n_samples = np.random.normal(mean, sd, 10000)
    
    # Calculate shear stress for all samples at once
    tan_fi = math.tan(math.radians(fi))  # Convert fi to radians if it's in degrees
    numerator = 0.52 * tan_fi
    denominator = ((1 + (np.abs(n_samples) * A))**2) * (1 + (1/2.5) * ((np.abs(n_samples) * F)**2) * tan_fi)
    shear_stress = numerator / denominator
    
    # Add to our data lists: repeat the current Re 10,000 times
    all_Re.extend([Re] * 10000)
    all_shear_stress.extend(shear_stress)

# Convert lists to NumPy arrays for smoother plotting
all_Re = np.array(all_Re)
all_shear_stress = np.array(all_shear_stress)

# Step 3: Create the log-log scatter plot
plt.figure(figsize=(10, 6))
# Use low alpha and small point size to handle overplotting of 100k points
plt.scatter(all_Re, all_shear_stress, alpha=0.05, s=1)

# Set log scales for both axes
plt.xscale('log')
plt.yscale('log')

# Add labels and title for clarity
plt.xlabel('Reynolds Number (Re)')
plt.ylabel('Shear Stress')
plt.title('Reynolds Number vs Shear Stress (Log-Log Scale)')

# Add grid lines for better readability
plt.grid(True, which="both", ls="-", alpha=0.3)

plt.show()

Key Improvements Explained: #

• Vectorized Operations: Using NumPy arrays instead of nested loops cuts down computation time significantly for 100,000 points.
• Proper Data Alignment: We repeat each Re 10,000 times to ensure every shear stress value has a matching x-axis value.
• Log-Log Scale: plt.xscale('log') and plt.yscale('log') enforce the logarithmic axes you requested.
• Overplotting Fix: The alpha=0.05 and s=1 settings make overlapping points visible instead of a solid, uninformative blob.

Don't forget to replace the placeholder constants (mean, sd, fi, etc.) with your actual values—especially note that fi is converted from degrees to radians since math.tan() expects radians as input.

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