Производная кубического сплайна является квадратичным сплайном. В SciPy есть только встроенный метод для поиска корней кубического сплайна. Итак, есть два подхода:
- Используйте сплайн 4-й степени для интерполяции, чтобы можно было легко найти корни его производной.
- Используйте кубический сплайн (который часто предпочтительнее) и напишите собственную функцию для корней его производной.
Я описываю оба решения ниже.
сплайн 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=3)
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.
person
Community
schedule
16.05.2018
ax**3+bx**2+cx+d
имеет свои (локальные) экстремумы на-b+-sqrt(b**2-3ac)/3a
- person Ma0   schedule 16.05.2018f
для разных значений, увеличивающихся наnp.finfo(float).eps
или несколько кратных им. - person Mohammad Athar   schedule 16.05.2018