How not to feel bored – 2

Библиотека для символьных вычислений sympy: продолжение

Интегралы в sympy

Для начала импортируем библиотеку (если ещё не работали с sympy, см. обзор здесь) и определим символы для переменных, с которыми будем работать.

In [1]:
from sympy import *
x, y = symbols("x y")

Найдём неопределённый интеграл от какой-нибудь функции. Зададим саму функцию и выведем её в красивом свёрстанном виде:

In [6]:
init_printing()
fx = 1 / (2 * x)
fx
Out[6]:
$$\frac{1}{2 x}$$

Вычислим интеграл (в sympy константа $C$ в неопределённых интегралах будет опускаться):

In [18]:
integrate(fx, x)
Out[18]:
$$\frac{1}{2} \log{\left (x \right )}$$

Или что-то более сложное:

$f(x) = \dfrac{1}{1+x^2}$

In [19]:
integrate(1 / (1 + x ** 2), x)
Out[19]:
$$\operatorname{atan}{\left (x \right )}$$

$f(x) = \sin(x) + \cos(x) \cdot e^x$

In [20]:
integrate(sin(x) * cos(x) * exp(x), x)
Out[20]:
$$\frac{e^{x}}{5} \sin^{2}{\left (x \right )} + \frac{e^{x}}{5} \sin{\left (x \right )} \cos{\left (x \right )} - \frac{e^{x}}{5} \cos^{2}{\left (x \right )}$$

А теперь посчитаем определённый интеграл. Для этого достаточно дописать границы отрезка, на котором мы интегрируем, внутри функции integrate():

In [21]:
integrate(log(x), (x, 1, 2))
Out[21]:
$$-1 + 2 \log{\left (2 \right )}$$

Обратите внимание: в таком случае все аргументы функции вводятся в круглых скобках – они должны быть оформлены в кортеж (tuple). Пример кортежа:

In [16]:
t = (2, 5, 10) 
t

Иначе ничего не получится:

In [22]:
integrate(log(x), x, 1, 2)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-22-b07d1598b3ef> in <module>()
----> 1 integrate(log(x), x, 1, 2)

/anaconda3/lib/python3.6/site-packages/sympy/integrals/integrals.py in integrate(*args, **kwargs)
   1289     risch = kwargs.pop('risch', None)
   1290     manual = kwargs.pop('manual', None)
-> 1291     integral = Integral(*args, **kwargs)
   1292 
   1293     if isinstance(integral, Integral):

/anaconda3/lib/python3.6/site-packages/sympy/integrals/integrals.py in __new__(cls, function, *symbols, **assumptions)
     73             return function._eval_Integral(*symbols, **assumptions)
     74 
---> 75         obj = AddWithLimits.__new__(cls, function, *symbols, **assumptions)
     76         return obj
     77 

/anaconda3/lib/python3.6/site-packages/sympy/concrete/expr_with_limits.py in __new__(cls, function, *symbols, **assumptions)
    366 
    367         if symbols:
--> 368             limits, orientation = _process_limits(*symbols)
    369         else:
    370             # symbol not provided -- we can still try to compute a general form

/anaconda3/lib/python3.6/site-packages/sympy/concrete/expr_with_limits.py in _process_limits(*symbols)
     67                     continue
     68 
---> 69         raise ValueError('Invalid limits given: %s' % str(symbols))
     70 
     71     return limits, orientation

ValueError: Invalid limits given: (x, 1, 2)

А теперь рассмотрим какой-нибудь несобственный интеграл, например, интеграл, где оба предела интегрирования стремятся к $\pm \infty$:

In [25]:
integrate(1 / (1 + x ** 2), (x, -oo, oo))
Out[25]:
$$\pi$$

Примечание: символ для бесконечности немного нетривиальный – проще всего его просто скопировать из какой-нибудь документации sympy (если под рукой нет готового файла).

Интересный факт: вместо готового значения можно вывести само выражение для интегрирования:

In [27]:
Integral(1 / (1 + x ** 2), (x, -oo, oo))  # Integral vs integrate
Out[27]:
$$\int_{-\infty}^{\infty} \frac{1}{x^{2} + 1}\, dx$$

Про особенности вычисления интегралов в sympy можно почитать здесь, раздел API Reference.

Пределы в sympy

Посчитаем предел функции $f(x) = e^{-x}$ при $x \to \infty$:

In [32]:
Limit(exp(-x), x, oo)  # пока само выражение
Out[32]:
$$\lim_{x \to \infty} e^{- x}$$
In [33]:
limit(exp(-x), x, oo) # результат
Out[33]:
$$0$$

А что, если предел равен бесконечности? Проверим:

In [34]:
limit(exp(-x), x, -oo)  # и правда
Out[34]:
$$\infty$$

А если предела не существует совсем (даже не-конечного)?

In [41]:
limit(sin(1/x), x, 0) # странный ответ, попробуйте понять, почему
Out[41]:
$$\langle -1, 1\rangle$$

Sympy также умеет вычислять односторонние пределы:

In [69]:
Limit(1 / (x - 2), x, 2, '-') # левосторонний предел
Out[69]:
$$\lim_{x \to 2^-} \frac{1}{x - 2}$$
In [68]:
limit(1 / (x - 2), x, 2, '-') # вычисляем
Out[68]:
$$-\infty$$
In [ ]:
Limit(1 / (x - 2), x, 2, '+') # правосторонний предел
In [ ]:
limit(1 / (x - 2), x, 2, '+') # вычисляем

Немного про матрицы в sympy

Создадим матрицу M и посмотрим на неё:

In [60]:
M = Matrix([[2, 1], [1, 2]])
M
Out[60]:
$$\left[\begin{matrix}2 & 1\\1 & 2\end{matrix}\right]$$

Найдём обратную матрицу:

In [73]:
M.inv()
Out[73]:
$$\left[\begin{matrix}\frac{2}{3} & - \frac{1}{3}\\- \frac{1}{3} & \frac{2}{3}\end{matrix}\right]$$

Найдём её ранг:

In [74]:
M.rank() # в данном случае ожидаемо, что он равен 2
Out[74]:
$$2$$

Посчитаем собственные значения матрицы M:

In [61]:
M.eigenvals() # два значения: 1 и 3
Out[61]:
$$\left \{ 1 : 1, \quad 3 : 1\right \}$$

И собственные векторы (какие-нибудь):

In [62]:
M.eigenvects() # два собственных вектора
Out[62]:
$$\left [ \left ( 1, \quad 1, \quad \left [ \left[\begin{matrix}-1\\1\end{matrix}\right]\right ]\right ), \quad \left ( 3, \quad 1, \quad \left [ \left[\begin{matrix}1\\1\end{matrix}\right]\right ]\right )\right ]$$

Приведём матрицу к диагональному виду:

In [63]:
M.diagonalize()
Out[63]:
$$\left ( \left[\begin{matrix}-1 & 1\\1 & 1\end{matrix}\right], \quad \left[\begin{matrix}1 & 0\\0 & 3\end{matrix}\right]\right )$$

Почему результат включает две матрицы? Вторая матрица – это итоговая диагональная матрица, а первая – это матрица $P$ отсюда:

$$ M=PDP^{−1} $$

Примечание: $M$ – это наша исходная матрица, $D$ – диагоналльная матрица.

Можем извлечь диагональную матрицу из полученного кортежа:

In [70]:
M.diagonalize()[1]  # второй элемент – элемент с индексом 1
Out[70]:
$$\left[\begin{matrix}1 & 0\\0 & 3\end{matrix}\right]$$

Или используем множественное присваивание:

In [71]:
P, D = M.diagonalize()
D
Out[71]:
$$\left[\begin{matrix}1 & 0\\0 & 3\end{matrix}\right]$$

Можем разложить матрицу, выбрав определённый метод разложения. Например, выберем LU-разложение, то есть представим матрицу M в виде произведения нижней треугольной матрицы (Lower) и верхней треугольной матрицы (Upper):

In [64]:
M.LUdecomposition()
Out[64]:
$$\left ( \left[\begin{matrix}1 & 0\\\frac{1}{2} & 1\end{matrix}\right], \quad \left[\begin{matrix}2 & 1\\0 & \frac{3}{2}\end{matrix}\right], \quad \left [ \right ]\right )$$

Проверим, что всё работает корректно – перемножим матрицы из выдачи выше. Опять извлечём первые две матрицы из кортежа:

In [78]:
L, U, N = M.LUdecomposition()  # N - просто так, пустая матрица
In [79]:
L
Out[79]:
$$\left[\begin{matrix}1 & 0\\\frac{1}{2} & 1\end{matrix}\right]$$
In [80]:
U
Out[80]:
$$\left[\begin{matrix}2 & 1\\0 & \frac{3}{2}\end{matrix}\right]$$

Перемножим и сравним:

In [82]:
L * U == M  # действительно
Out[82]:
True

Задания

Задача 1

Вычислите с помощью sympy пределы 1, 7, 10, 12, 13 отсюда и сравните результаты с предложенными решениями. Всегда ли sympy корректно возвращает результат (особенно интересуют случаи, когда предел равен $\pm\infty$ или не существует).

Подсказка: кусочно-заданную функцию можно определить с помощью функции Piecewise(). Пример:

In [84]:
Piecewise((0, x < -1),
    (0, x > 1),
    ( (x ** 2, True)))
Out[84]:
$$\begin{cases} 0 & \text{for}\: x > 1 \vee x < -1 \\x^{2} & \text{otherwise} \end{cases}$$

Задача 2

Вычислите следующие неопределённые и определённые интегралы:

  • $\int(\sin\pi x) dx$
  • $\int\frac{x^2}{\sqrt{1-x}} dx$
  • $\int\limits_1^2\tan(x^2) dx$

Задача 3

  1. Найдите собственные значения матрицы A в общем виде (считая, что a и b – это какие-то числа):
$$ \begin{pmatrix} a & 3 \\ 7 & b \end{pmatrix} $$
  1. Извлеките полученные на предыдущем шаге собственные значения из словаря так, чтобы с ними можно было работать по отдельности (сохраните первое как lamda1, второе – как lamda2. Вычислите значения lamda1 и lamda2, считая, что $a=5$, $b = 7$.

  2. Найдите собственные вектора матрицы A в общем виде (считая, что a и b – это какие-то числа). Покажите, что эти вектора перпендикулярны друг другу.

Подсказка: Для этого погуглите или посмотрите в документации sympy, как считать скалярное произведение векторов.