Использование удобных классов

Удобные классы, предоставляемые пакетом полиномов, это:

Ряды в данном контексте представляют собой конечные суммы соответствующих полиномиальных базисных функций, умноженных на коэффициенты. Например, степенной ряд выглядит следующим образом:

и имеет коэффициенты . Ряд Чебышева с теми же коэффициентами выглядит следующим образом:

и более общим образом

где в данном случае — это функции Чебышева степени , но они могут быть и базисными функциями любого другого класса. Согласно соглашению для всех классов, коэффициент соответствует базисной функции степени i.

Все классы являются неизменяемыми и имеют одинаковые методы, а именно, они реализуют числовые операторы Python +, -, *, //, %, divmod, **, == и !=. Два последних оператора могут быть немного проблематичными из-за ошибок округления с плавающей точкой. Теперь мы представим краткую демонстрацию различных операций, используя NumPy версии 1.7.0.

Основы

Сначала нам нужен класс полиномов и экземпляр полинома для работы. Классы можно импортировать непосредственно из пакета полиномов или из модуля соответствующего типа. Здесь мы импортируем из пакета и используем стандартный класс Polynomial из-за его знакомой функциональности:

>>> from numpy.polynomial import Polynomial as P
>>> p = P([1,2,3])
>>> p
Polynomial([ 1.,  2.,  3.], domain=[-1,  1], window=[-1,  1])

Обратите внимание, что в длинной версии вывода есть три части. Первая — коэффициенты, вторая — область определения, а третья — интервал:

>>> p.coef
array([ 1.,  2.,  3.])
>>> p.domain
array([-1.,  1.])
>>> p.window
array([-1.,  1.])

Вывод полинома имеет более короткий вид без области определения и интервала:

>>> print p
poly([ 1.  2.  3.])

Мы рассмотрим область определения и интервал, когда дойдем до подгонки, а пока мы их игнорируем и перейдём к основным алгебраическим и арифметическим операциям.

Сложение и вычитание:

>>> p + p
Polynomial([ 2.,  4.,  6.], domain=[-1,  1], window=[-1,  1])
>>> p - p
Polynomial([ 0.], domain=[-1,  1], window=[-1,  1])

Умножение:

>>> p * p
Polynomial([  1.,   4.,  10.,  12.,   9.], domain=[-1,  1], window=[-1,  1])

Степень:

>>> p**2
Polynomial([  1.,   4.,  10.,  12.,   9.], domain=[-1,  1], window=[-1,  1])

Деление:

Целочисленное деление, ‘//’, — это оператор деления для классов полиномов. Полиномы в этом отношении обрабатываются как целые числа. Для версий Python < 3.x оператор ‘/’ отображается на ‘//’, как и в Python, а для более поздних версий оператор ‘/’ будет работать только для деления на скаляры. В какой-то момент он будет устаревшим:

>>> p // P([-1, 1])
Polynomial([ 5.,  3.], domain=[-1,  1], window=[-1,  1])

Остаток:

>>> p % P([-1, 1])
Polynomial([ 6.], domain=[-1,  1], window=[-1,  1])

Divmod:

>>> quo, rem = divmod(p, P([-1, 1]))
>>> quo
Polynomial([ 5.,  3.], domain=[-1,  1], window=[-1,  1])
>>> rem
Polynomial([ 6.], domain=[-1,  1], window=[-1,  1])

Вычисление:

>>> x = np.arange(5)
>>> p(x)
array([  1.,   6.,  17.,  34.,  57.])
>>> x = np.arange(6).reshape(3,2)
>>> p(x)
array([[  1.,   6.],
       [ 17.,  34.],
       [ 57.,  86.]])

Подстановка:

Подставьте полином для x и разложите результат. Здесь мы подставляем p в себя, что приводит к новому полиному степени 4 после разложения. Если полиномы рассматривать как функции, это композиция функций:

>>> p(p)
Polynomial([  6.,  16.,  36.,  36.,  27.], domain=[-1,  1], window=[-1,  1])

Корни:

>>> p.roots()
array([-0.33333333-0.47140452j, -0.33333333+0.47140452j])

Не всегда удобно использовать экземпляры Polynomial явно, поэтому кортежи, списки, массивы и скаляры автоматически преобразуются в арифметических операциях:

>>> p + [1, 2, 3]
Polynomial([ 2.,  4.,  6.], domain=[-1,  1], window=[-1,  1])
>>> [1, 2, 3] * p
Polynomial([  1.,   4.,  10.,  12.,   9.], domain=[-1,  1], window=[-1,  1])
>>> p / 2
Polynomial([ 0.5,  1. ,  1.5], domain=[-1,  1], window=[-1,  1])

Полиномы, отличающиеся по области определения, интервалу или классу, не могут быть смешаны в арифметических операциях:

>>> from numpy.polynomial import Chebyshev as T
>>> p + P([1], domain=[0,1])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<string>", line 213, in __add__
TypeError: Domains differ
>>> p + P([1], window=[0,1])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<string>", line 215, in __add__
TypeError: Windows differ
>>> p + T([1])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<string>", line 211, in __add__
TypeError: Polynomial types differ

Но разные типы могут быть использованы для подстановки. Фактически, именно так происходит преобразование классов Polynomial друг в друга для преобразования типа, области определения и интервала:

>>> p(T([0, 1]))
Chebyshev([ 2.5,  2. ,  1.5], domain=[-1,  1], window=[-1,  1])

Что даёт полином p в форме Чебышева. Это работает, потому что и подстановка на не изменяет исходный полином. Однако все умножения и деления будут выполняться с использованием рядов Чебышева, поэтому тип результата будет таким.

Предполагается, что все экземпляры полиномов являются неизменяемыми, поэтому расширенные операции (+=, -=, и т.д.) и любая другая функциональность, которая нарушала бы неизменяемость экземпляра полинома, намеренно не реализована.

Дифференциальное исчисление

Экземпляры полиномов могут быть интегрированы и дифференцированы:

>>> from numpy.polynomial import Polynomial as P
>>> p = P([2, 6])
>>> p.integ()
Polynomial([ 0.,  2.,  3.], domain=[-1,  1], window=[-1,  1])
>>> p.integ(2)
Polynomial([ 0.,  0.,  1.,  1.], domain=[-1,  1], window=[-1,  1])

В первом примере p интегрируется один раз, а во втором примере интегрируется дважды. По умолчанию нижняя граница интегрирования и постоянная интегрирования равны 0, но оба значения могут быть указаны:

>>> p.integ(lbnd=-1)
Polynomial([-1.,  2.,  3.], domain=[-1,  1], window=[-1,  1])
>>> p.integ(lbnd=-1, k=1)
Polynomial([ 0.,  2.,  3.], domain=[-1,  1], window=[-1,  1])

В первом случае нижняя граница интегрирования установлена в -1, а постоянная интегрирования равна 0. Во втором случае постоянная интегрирования также установлена в 1. Дифференцирование проще, так как единственный параметр — это количество раз, которое полином дифференцируется:

>>> p = P([1, 2, 3])
>>> p.deriv(1)
Polynomial([ 2.,  6.], domain=[-1,  1], window=[-1,  1])
>>> p.deriv(2)
Polynomial([ 6.], domain=[-1,  1], window=[-1,  1])

Другие конструкторы полиномов

Создание полиномов, задавая коэффициенты, — это лишь один способ получения экземпляра полинома. Они также могут быть созданы, задавая их корни, путем преобразования из других типов полиномов и методом наименьших квадратов. Подгонка обсуждается в отдельном разделе, а другие методы продемонстрированы ниже:

>>> from numpy.polynomial import Polynomial as P
>>> from numpy.polynomial import Chebyshev as T
>>> p = P.fromroots([1, 2, 3])
>>> p
Polynomial([ -6.,  11.,  -6.,   1.], domain=[-1,  1], window=[-1,  1])
>>> p.convert(kind=T)
Chebyshev([ -9.  ,  11.75,  -3.  ,   0.25], domain=[-1,  1], window=[-1,  1])

Метод convert также может преобразовывать область определения и интервал:

>>> p.convert(kind=T, domain=[0, 1])
Chebyshev([-2.4375 ,  2.96875, -0.5625 ,  0.03125], [ 0.,  1.], [-1.,  1.])
>>> p.convert(kind=P, domain=[0, 1])
Polynomial([-1.875,  2.875, -1.125,  0.125], [ 0.,  1.], [-1.,  1.])

В версиях numpy >= 1.7.0 также доступны классы методов basis и cast. Метод cast работает как метод convert, а метод basis возвращает полином базиса заданной степени:

>>> P.basis(3)
Polynomial([ 0.,  0.,  0.,  1.], domain=[-1,  1], window=[-1,  1])
>>> T.cast(p)
Chebyshev([ -9.  ,  11.75,  -3.  ,   0.25], domain=[-1,  1], window=[-1,  1])

Преобразования между типами могут быть полезны, но их не рекомендуется использовать для рутинной работы. Потеря точности чисел при переходе от ряда Чебышева степени 50 к ряду полиномов той же степени может сделать результаты числового вычисления по существу случайными.

Подгонка

Подгонка — это причина, по которой атрибуты domain и window являются частью удобных классов. Для иллюстрации проблемы ниже показаны значения полиномов Чебышева до степени 5.

>>> import matplotlib.pyplot as plt
>>> from numpy.polynomial import Chebyshev as T
>>> x = np.linspace(-1, 1, 100)
>>> for i in range(6): ax = plt.plot(x, T.basis(i)(x), lw=2, label="$T_%d$"%i)
...
>>> plt.legend(loc="upper left")
<matplotlib.legend.Legend object at 0x3b3ee10>
>>> plt.show()

В диапазоне -1 <= x <= 1 они представляют собой хорошие, равнопеременные функции, лежащие между +/- 1. Те же графики в диапазоне -2 <= x <= 2 выглядят совершенно иначе:

>>> import matplotlib.pyplot as plt
>>> from numpy.polynomial import Chebyshev as T
>>> x = np.linspace(-2, 2, 100)
>>> for i in range(6): ax = plt.plot(x, T.basis(i)(x), lw=2, label="$T_%d$"%i)
...
>>> plt.legend(loc="lower right")
<matplotlib.legend.Legend object at 0x3b3ee10>
>>> plt.show()

Как видно, «хорошие» части уменьшились до незначительности. При использовании полиномов Чебышева для подгонки мы хотим использовать область, где x находится между -1 и 1, и именно это задает window. Однако маловероятно, что данные, которые подлежат подгонке, содержат все точки данных в этом интервале, поэтому мы используем domain для указания интервала, где находятся точки данных. При выполнении подгонки область определения сначала отображается на интервал с помощью линейного преобразования, и выполняется обычная подгонка методом наименьших квадратов, используя сопоставленные точки данных. Интервал и область определения подгонки являются частью возвращаемого ряда и автоматически используются при вычислении значений, производных и т.д. Если они не указаны в вызове, процедура подгонки будет использовать стандартный интервал и наименьшую область определения, которая содержит все точки данных. Это проиллюстрировано ниже для подгонки к шумовой синусоиде.

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from numpy.polynomial import Chebyshev as T
>>> np.random.seed(11)
>>> x = np.linspace(0, 2*np.pi, 20)
>>> y = np.sin(x) + np.random.normal(scale=.1, size=x.shape)
>>> p = T.fit(x, y, 5)
>>> plt.plot(x, y, 'o')
[<matplotlib.lines.Line2D object at 0x2136c10>]
>>> xx, yy = p.linspace()
>>> plt.plot(xx, yy, lw=2)
[<matplotlib.lines.Line2D object at 0x1cf2890>]
>>> p.domain
array([ 0.        ,  6.28318531])
>>> p.window
array([-1.,  1.])
>>> plt.show()