Discrete Integration with Overlapping Partitions of Fixed Width

# Discrete Integration with Overlapping Partitions of Fixed Width
# Uses:
# LP Positions on an AMM
# Solid to Liquid transition Ice >> Water
# Bose-Einstein condensate converting to Fermi-Dirac statistics

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf

def f(x):
    return np.exp(-((x - 10)**2) / 10)

a, b = 3, 19  # Discrete Integration Bounds
width = 3 # our constraint
step = 0.5 # our discrete grid
x_starts = np.arange(a-step, b + step, step)

# Build LP Positions
rectangles = []
for x in x_starts:
    base = sum(r[2] for r in rectangles if r[0] <= x < r[1])
    height = max(0, f(x) - base)
    rectangles.append((x, x + width, height, base))

# Compute total area
total_area = sum(height * width for _, _, height, _ in rectangles)
exact_area = np.sqrt(10 * np.pi) / 2 * (erf(9 / np.sqrt(10)) + erf(7 / np.sqrt(10)))

# Cumulative height data
x_points = np.arange(min(x_starts), max(x_starts) + width, step)
heights = [sum(r[2] for r in rectangles if r[0] <= x < r[1]) for x in x_points]

# Midpoint Shift for plotting only
shift = step / 2
x_points_shifted = x_points - shift

plt.figure(figsize=(14, 7))
x_curve = np.linspace(a - 3, b + 3, 200)
y_curve = f(x_curve)
plt.plot(x_curve, y_curve, 'b-', label='f(x) = e^(-(x-10)²/10)', linewidth=2.5)
plt.step(x_points_shifted, heights, where='post', color='r', label='Cumulative Height (Shifted Step/2)', linewidth=2)

# Rectangles from y=0 (shifted)
for i, (x0, x1, height, _) in enumerate(rectangles):
    if height > 0:
        color = plt.cm.viridis(i / len(rectangles))
        plt.gca().add_patch(
            plt.Rectangle((x0 - shift, 0), width, height, alpha=0.6, facecolor=color,
                          edgecolor='k', linewidth=0.8,
                          label=f'LP Positions' if i == 0 else None)
        )

plt.title(f'Centered Stacking with Midpoint\nApprox Area = {total_area:.4f}, Exact = {exact_area:.4f}',
          fontsize=14, pad=10)
plt.xlabel('x', fontsize=12)
plt.ylabel('y', fontsize=12)
plt.xlim(a - 3, b + 3)
plt.ylim(0, 1.2)
plt.grid(True, linestyle='--', alpha=0.5)
plt.legend(loc='upper left', fontsize=10)
plt.tight_layout()
plt.show()

print("Sample rectangles (unshifted positions):")
for x0, x1, h, base in rectangles[:5]:
    print(f"Rect from {x0} to {x1}: height = {h:.4f}, base = {base:.4f}, top = {base + h:.4f}")
print(f"Rectangles near peak (unshifted positions):")
for x0, x1, h, base in rectangles[15:20]:
    print(f"Rect from {x0} to {x1}: height = {h:.4f}, base = {base:.4f}, top = {base + h:.4f}")
print(f"Total area: {total_area:.4f}, Exact: {exact_area:.4f}")