User:Sintel/Zeta plot python

From Xenharmonic Wiki
Revision as of 20:41, 8 April 2025 by Sintel (talk | contribs) (Created page with "Code for plotting some zetas <syntaxhighlight lang="python"> import math import numpy as np import matplotlib.pyplot as plt from mpmath import zeta, mp from scipy.signal import find_peaks # Set precision for mpmath mp.dps = 15 # Define the range for t t_min, t_max = 4.8, 42 # Real part of the variable # Use 0.5 for critical line sigma = 0.5 # Create array of t values num_points = 2000 t_values = np.linspace(t_min, t_max, num_points) # Calculate zeta function magni...")
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

Code for plotting some zetas


import math
import numpy as np
import matplotlib.pyplot as plt
from mpmath import zeta, mp
from scipy.signal import find_peaks

# Set precision for mpmath
mp.dps = 15

# Define the range for t
t_min, t_max = 4.8, 42

# Real part of the variable
# Use 0.5 for critical line
sigma = 0.5

# Create array of t values
num_points = 2000
t_values = np.linspace(t_min, t_max, num_points)

# Calculate zeta function magnitude using s = sigma + 2πt/ln(2)
zeta_magnitudes = np.zeros(num_points)
for i, t in enumerate(t_values):
    s = complex(sigma, 2 * math.pi * t / math.log(2))
    zeta_magnitudes[i] = float(abs(zeta(s)))

# Create the plot
plt.figure(figsize=(24, 7))
plt.plot(t_values, zeta_magnitudes, color="black")

# Add title and labels
plt.title(f"Magnitude of Riemann Zeta Function for s = {sigma:.2f} + 2πt/ln(2)", fontsize=14)
plt.xlabel("t", fontsize=12)
plt.ylabel("|ζ(s)|", fontsize=12)

# Add grid lines
# plt.grid(alpha=0.3)

# Find local maxima
peaks, _ = find_peaks(zeta_magnitudes, height=0.5, distance=10)
peak_t_values = t_values[peaks]
peak_magnitudes = zeta_magnitudes[peaks]

# Find record maxima (peaks that are higher than any previous peak)
record_peaks = []
record_t_values = []
record_magnitudes = []
max_magnitude_so_far = 0

for i, (t, mag) in enumerate(zip(peak_t_values, peak_magnitudes)):
    if mag > max_magnitude_so_far:
        max_magnitude_so_far = mag
        record_peaks.append(i)
        record_t_values.append(t)
        record_magnitudes.append(mag)

# Mark all record peaks
plt.plot(record_t_values, record_magnitudes, "ro", markersize=4)

for t, mag in zip(record_t_values, record_magnitudes):
    # Add text label above
    name = f"{round(t)}edo"
    plt.text(
        t,
        mag + 0.25,
        name,
        ha="center",
        fontsize=10,
    )

plt.axhline(y=0, color="black", alpha=0.2)

plt.tight_layout()

ax = plt.gca()
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)

plt.subplots_adjust(bottom=0.15)  # Make room for the explanation text
plt.show()