Finding the maximum value of a curve

Question

I have two lists, which I have plotted on a graph:

y = [-0.44375618032407804, -0.4330176209985608, -0.4221413824497001, -0.4111868108222988, -0.4002148653666401, -0.3892876774098184, -0.3784680960581614, -0.36781922437211045, -0.3574039498559348, -0.34728447317173006, -0.33752183901958394, -0.32817547312320006, -0.3193027292225641, -0.31095844990269483, -0.30319454498070525, -0.29605959103321394, -0.28959845547379326, -0.28385194838712835, -0.2788565050946762, -0.27464390216791346, -0.2712410093220698, -0.26866957931806573, -0.2669460776759942, -0.26608155366278724, -0.26608155366278724, -0.2669460776759942, -0.26866957931806573, -0.2712410093220698, -0.27464390216791346, -0.2788565050946762, -0.28385194838712835, -0.28959845547379326, -0.29605959103321394, -0.30319454498070525, -0.31095844990269483, -0.3193027292225641, -0.32817547312320006, -0.33752183901958394, -0.34728447317173006, -0.3574039498559348, -0.36781922437211045, -0.37846809605816145, -0.3892876774098184, -0.4002148653666401, -0.41118681082229874, -0.4221413824497001, -0.4330176209985608, -0.44375618032407793]

x = [0.10285714285714286, 0.15428571428571428, 0.2057142857142857, 0.2571428571428571, 0.30857142857142855, 0.36, 0.4114285714285714, 0.46285714285714286, 0.5142857142857142, 0.5657142857142857, 0.6171428571428571, 0.6685714285714286, 0.72, 0.7714285714285715, 0.8228571428571428, 0.8742857142857142, 0.9257142857142857, 0.9771428571428571, 1.0285714285714285, 1.08, 1.1314285714285715, 1.1828571428571428, 1.2342857142857142, 1.2857142857142856, 1.3371428571428572, 1.3885714285714286, 1.44, 1.4914285714285713, 1.542857142857143, 1.5942857142857143, 1.6457142857142857, 1.697142857142857, 1.7485714285714284, 1.8, 1.8514285714285714, 1.9028571428571428, 1.9542857142857142, 2.005714285714286, 2.057142857142857, 2.1085714285714285, 2.16, 2.2114285714285713, 2.262857142857143, 2.314285714285714, 2.3657142857142857, 2.4171428571428573, 2.4685714285714284, 2.52]

I am looking to find the coordinates of the highest point of the curve produced (not necessarily the max values of x and y). I can see from other similar questions that scipy's curve_fit and minimize_scalar functions are effective in doing this, but I cannot understand how to properly apply them when I already have the required data points. It looks like one of the parameters for curve_fit is a function that transforms the data in some way - I do not want to do this.

I am likely misunderstanding something fundamental about how this works...

Answer 1

The following approach fits a circle through the given point with the highest y-value and the two surrounding points. Then it takes the highest point of that circle (being the x of the center and its y plus the radius).

The formulas were calculated via sympy.

The plot shows the difference between just taking the highest point and taking the highest point of the fitted circle. The y-values are quite close, the x-values differ. If nothing is known about the underlying curve, fitting a circle would be closer to reality than just taking the maximum, especially concerning the x-value. Of course, the real maximum can't be known.

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

def find_highest(x, y):
    ind = np.argmax(y)
    if ind == 0 or ind == len(y) - 1:
        return x[ind], y[ind]
    else:
        x0, x1, x2 = x[ind - 1], x[ind], x[ind + 1]
        y0, y1, y2 = y[ind - 1], y[ind], y[ind + 1]
        d = x0 * y1 - x0 * y2 - x1 * y0 + x1 * y2 + x2 * y0 - x2 * y1
        if abs(d) < 1e-9:
            return x[ind], y[ind]
        else:
            e = x0 * x1 - x0 * x2 - x1 * x2 + x2 ** 2 + y0 * y1 - y0 * y2 - y1 * y2 + y2 ** 2
            x_center = ((x0 + x1) + (y0 - y1) * e / d) / 2
            y_center = ((x1 - x0) * e / d + (y0 + y1)) / 2
            rad = sqrt(
                (x0 - ((x0 + x1) + (y0 - y1) * e / d) / 2) ** 2 + (y0 - ((x1 - x0) * e / d + (y0 + y1)) / 2) ** 2)
            return x_center, y_center + rad

# x = [0.10285714285714286, ...]
# y = [-0.44375618032407804, ...]

plt.plot(x, y)
ind = np.argmax(y)
print(f"Highest point using only the given points: ({x[ind]}, {y[ind]})")
plt.plot(x[ind], y[ind], 'g.')

xh, yh = find_highest(x, y)
plt.plot(xh, yh, 'ro')
plt.axhline(yh, color='r', ls=':')
print(f"Highest point using fitted circle: ({xh}, {yh})")
plt.show()

Highest point using only the given points: (1.2857142857142856, -0.26608155366278724)
Highest point using fitted circle:         (1.3114285714285683, -0.26597350533539954)

Zoomed-in plot:

The sympy code to calculate the center and the radius:

from sympy import symbols
from sympy.geometry import Point, Circle
from sympy.abc import d, e

x0, y0, x1, y1, x2, y2 = symbols('x0, y0, x1, y1, x2, y2', real=True)
C0 = Circle(Point(x0, y0), Point(x1, y1), Point(x2, y2))
subs_dict = {x0 * y1 - x0 * y2 - x1 * y0 + x1 * y2 + x2 * y0 - x2 * y1: d,
             x0 * x1 - x0 * x2 - x1 * x2 + x2 ** 2 + y0 * y1 - y0 * y2 - y1 * y2 + y2 ** 2: e}
print('center:', C0.center)
print('radius:', C0.radius)

Here two common terms can be separated:

subs_dict = {x0 * y1 - x0 * y2 - x1 * y0 + x1 * y2 + x2 * y0 - x2 * y1: d,
             x0 * x1 - x0 * x2 - x1 * x2 + x2 ** 2 + y0 * y1 - y0 * y2 - y1 * y2 + y2 ** 2: e}
print('center:', C0.center.subs(subs_dict))
print('radius:', C0.radius.subs(subs_dict))

A similar approach could fit a quadratic curve through the highest y-value and the two surrounding points. Then, the highest point on the quadratic curve can be found by setting the derivative to zero. This returns the same highest point (up till 6 decimals).

def find_highest_via_quadratic(x, y):
    ind = np.argmax(y)
    if ind == 0 or ind == len(y) - 1:
        return x[ind], y[ind]
    else:
        x0, x1, x2 = x[ind - 1], x[ind], x[ind + 1]
        y0, y1, y2 = y[ind - 1], y[ind], y[ind + 1]
        d = (x0 - x1) * (x0 - x2) * (x1 - x2)
        a = (x0 * (y2 - y1) + x1 * (y0 - y2) + x2 * (y1 - y0))
        b = (x0 ** 2 * (y1 - y2) + x1 ** 2 * (y2 - y0) + x2 ** 2 * (y0 - y1))
        c = (x0 ** 2 * (x1 * y2 - x2 * y1) + x1 ** 2 * (x2 * y0 - x0 * y2) + x2 ** 2 * (x0 * y1 - x1 * y0))
        if abs(d) < 1e-9 or abs(a) < 1e-9:
            return x[ind], y[ind]
        else:
            x_high = -b / (2 * a)
            y_high = ((a * x_high + b) * x_high + c) / d
            return x_high, y_high

Answer 2

I don't really understand why this wouldn't be the max value of the array y.

So all you'd do is:

max_id = np.argmax(y)
print(f"The Peak is at (x,y)=({x[max_id]},{y[max_id]})")