Нахождение максимального значения кривой

Question

У меня есть два списка, которые я построил на графике:

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]

введите описание изображения здесь

Я ищу координаты наивысшей точки полученной кривой (не обязательно максимальные значения x и y). Из других подобных вопросов я вижу, что функции scipy curve_fit и minimize_scalar эффективны при этом, но я не могу понять, как правильно их применять, когда у меня уже есть необходимые точки данных. Похоже, что один из параметров для curve_fit - это функция, которая каким-то образом преобразует данные - я не хочу этого делать.

Я, вероятно, неправильно понимаю что-то фундаментальное о том, как это работает ...

Answer 1

Следующий подход помещает окружность через заданную точку с наивысшим значением y и две окружающие точки. Затем он берет наивысшую точку этого круга (являющуюся x центра и его y плюс радиус).

Формулы были рассчитаны с помощью sympy.

График показывает разницу между взятие наивысшей точки и взятие наивысшей точки подобранного круга. Значения y довольно близки, значения x отличаются. Если о лежащей в основе кривой ничего не известно, построение круга было бы ближе к реальности, чем просто взятие максимума, особенно в отношении значения x. Конечно, реальный максимум неизвестен.

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)

Увеличенный график:

The sympy код для вычисления центра и радиус:

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)

Здесь можно разделить два общих термина:

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))

Аналогичный подход может соответствовать квадратичной кривой c через наивысшее значение y и два окружающие точки. Затем можно найти наивысшую точку на кривой квадратичной c, установив производную равной нулю. Это возвращает ту же самую высокую точку (до 6 знаков после запятой).

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

Я действительно не понимаю, почему это не максимальное значение массива y.

Итак, все, что вам нужно сделать, это:

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