User:Sintel/Zeta plot python

From Xenharmonic Wiki
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πit/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πit/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.show()