Найти максимум / минимум 1d интерполированной функции - PullRequest
Новая
       51

Найти максимум / минимум 1d интерполированной функции

0 голосов
/ 16 мая 2018

У меня есть набор данных, которые я интерполирую с kind = 'cubic'.

Я хотел бы найти максимум этой кубической функции интерполяции.

В настоящее время я просто нахожу максимальное значение в массиве интерполированных данных, но мне было интересно, можно ли дифференцировать интерполированную функцию как объект, чтобы найти ее экстремумы?

Код: 

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

x_axis = np.array([ 2.14414414,  2.15270826,  2.16127238,  2.1698365 ,  2.17840062, 2.18696474,  2.19552886,  2.20409298,  2.2126571 ,  2.22122122])
y_axis = np.array([ 0.67958442,  0.89628424,  0.78904004,  3.93404167,  6.46422317, 6.40459954,  3.80216674,  0.69641825,  0.89675386,  0.64274198])

f = interp1d(x_axis, y_axis, kind = 'cubic')

x_new = np.linspace(x_axis[0], x_axis[-1],100)

fig = plt.subplots()
plt.plot(x_new, f(x_new))

1 Ответ

0 голосов
/ 16 мая 2018

Производная кубического сплайна является квадратичным сплайном. В SciPy есть только встроенный метод для поиска корней кубического сплайна. Итак, есть два подхода:

  1. Используйте сплайн 4-й степени для интерполяции, чтобы можно было легко найти корни его производной.
  2. Используйте кубический сплайн (который часто предпочтителен) и напишите пользовательскую функцию для корней его производной.

Я опишу оба решения ниже.

Сплайн 4-й степени

Использовать InterpolatedUnivariateSpline ., У которого есть .derivative метод, возвращающий кубический сплайн, к которому можно применить метод .roots. 

from scipy.interpolate import InterpolatedUnivariateSpline
f = InterpolatedUnivariateSpline(x_axis, y_axis, k=4)
cr_pts = f.derivative().roots()
cr_pts = np.append(cr_pts, (x_axis[0], x_axis[-1]))  # also check the endpoints of the interval
cr_vals = f(cr_pts)
min_index = np.argmin(cr_vals)
max_index = np.argmax(cr_vals)
print("Maximum value {} at {}\nMinimum value {} at {}".format(cr_vals[max_index], cr_pts[max_index], cr_vals[min_index], cr_pts[min_index]))

Выход:

Максимальное значение 6,779687224066201 при 2,1824928509277037
Минимальное значение 0.34588448400295346 при 2.2075868177297036

кубический сплайн

Нам нужна пользовательская функция для корней квадратичного сплайна. Вот оно (объяснено ниже). 

def quadratic_spline_roots(spl):
    roots = []
    knots = spl.get_knots()
    for a, b in zip(knots[:-1], knots[1:]):
        u, v, w = spl(a), spl((a+b)/2), spl(b)
        t = np.roots([u+w-2*v, w-u, 2*v])
        t = t[np.isreal(t) & (np.abs(t) <= 1)]
        roots.extend(t*(b-a)/2 + (b+a)/2)
    return np.array(roots)

Теперь действуйте точно так же, как описано выше, за исключением использования специального решателя. 

from scipy.interpolate import InterpolatedUnivariateSpline
f = InterpolatedUnivariateSpline(x_axis, y_axis, k=4)
cr_pts = quadratic_spline_roots(f.derivative())
cr_pts = np.append(cr_pts, (x_axis[0], x_axis[-1]))  # also check the endpoints of the interval
cr_vals = f(cr_pts)
min_index = np.argmin(cr_vals)
max_index = np.argmax(cr_vals)
print("Maximum value {} at {}\nMinimum value {} at {}".format(cr_vals[max_index], cr_pts[max_index], cr_vals[min_index], cr_pts[min_index]))

Выход:

Максимальное значение 6,782781181150518 при 2,1824928579767167
Минимальное значение 0,45017143148176136 при 2,2070746522580795

Небольшое расхождение с выводом в первом методе не является ошибкой; сплайн 4-й степени и сплайн 3-й степени немного отличаются.

Объяснение quadratic_spline_roots

Предположим, мы знаем, что значения квадратичного полинома в -1, 0, 1 равны u, v, w. Каковы его корни на интервале [-1, 1]? С некоторой алгеброй мы можем найти, что многочлен равен 

((u+w-2*v) * x**2 + (w-u) * x + 2*v) / 2

Теперь можно использовать квадратную формулу, но лучше использовать np.roots, поскольку она также будет обрабатывать случай, когда старший коэффициент равен нулю. Затем корни фильтруются до действительных чисел от -1 до 1. Наконец, если интервал равен [a, b] вместо [-1, 1], выполняется линейное преобразование.

Бонус: ширина кубического сплайна в среднем диапазоне

Предположим, мы хотим найти, где сплайн принимает значение, равное среднему значению его максимума и минимума (то есть его среднего диапазона). Тогда мы обязательно должны использовать кубический сплайн для интерполяции, потому что теперь для него понадобится метод roots. Нельзя просто сделать (f - mid_range).roots(), так как добавление константы в сплайн не поддерживается в SciPy. Вместо этого постройте смещенный вниз сплайн из y_axis - mid_range. 

mid_range = (cr_vals[max_index] + cr_vals[min_index])/2
f_shifted = InterpolatedUnivariateSpline(x_axis, y_axis - mid_range, k=3)
roots = f_shifted.roots()
print("Mid-range attained from {} to {}".format(roots.min(), roots.max()))

Средний уровень достигнут от 2,169076230034363 до 2,195974299834667

...