Найти максимум / минимум 1D интерполированной функции
У меня есть набор данных, которые я интерполирую с помощью 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 ответ:
Производная кубического сплайна является квадратичным сплайном. У SciPy есть только встроенный метод, чтобы найти корни кубического сплайна. Итак, есть два подхода:
- используйте сплайн 4-й степени для интерполяции, чтобы корни его производной можно было легко найти.
- используйте кубический сплайн (что часто предпочтительнее) и напишите пользовательскую функцию для корней его производной.
Я опишу оба решения ниже.
4-я степень сплайн
ИспользуйтеИнтерполяцию Univariatespline .который имел метод
.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]))
Вывод:
Небольшое расхождение с выводом в первом методе не является ошибкой; сплайн 4-й степени и сплайн 3-й степени немного отличаются.Максимальное значение 6.782781181150518 при 2.1824928579767167
Минимальное значение 0,45017143148176136 при 2,2070746522580795Объяснение
Предположим, мы знаем, что значения квадратичного полинома при -1, 0, 1 равны u, v, w. каковы его корни на интервале [-1, 1]? С помощью некоторой алгебры мы можем найти, что многочлен равенquadratic_spline_roots
Теперь можно использовать квадратичную формулу, но лучше использовать((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