Кратчайшее введение в Python для математики

Илья Щуров, НИУ ВШЭ.

Документ распространяется по лицензии CC BY-SA 4.0 Unported. Все фрагменты кода в этом блокноте переданы в общественное достояние. Исходные коды доступны на github.

Что это?

Этот документ предназначен для тех, кто хочет быстро начать использовать Python при решении математических и околоматематических задач. В нём предполагается, что вы уже обладаете базовыми навыками программирования, но не знаете Python. Если таких навыков нет, вы, возможно, сможете использовать отдельные рецепты отсюда, но вряд ли получите ту свободу, которую даёт знание универсального языка программирования.

Если вы никогда раньше не программировали — или наоборот, хотите освоить базовый Python более фундаментально, а уже потом переходить к его математическим приложениям — я рекомендую интерактивный курс Pythontutor.ru (там есть визуализатор пошагового выполнения кода и задачи с автоматической проверкой) или мой курс по Python в ВШЭ (там есть видеолекции, ко всему прочему).

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

Всё взять и запустить!

Самый простой способ запустить Python — зайти на http://tmpnb.org/ (в IE может не работать, лучше использовать Firefox или Chrome) и выбрать New → Python 3. Откроется notebook, состоящий из ячеек. Если в ячейку ввести код и нажать Shift + Enter, он выполнится.

Чуть сложнее: скачать Anaconda, установить, найти в меню «Пуск» или его аналоге Jupyter Notebook или IPython Notebook и запустить. (Под Mac OS это Anaconda Launcher.) Откроется окно браузера, такое же, как при открытии http://tmpnb.org/.

Что-что, простите?

Немного терминологии, чтобы не запутаться.

  • Python — это язык программирования. Он бывает двух основных модификаций — Python 2 и Python 3. Конкретная версия имеет обозначение через точку - например, текущая (на начало 2016 года) версия Python 3 имеет номер 3.5, а Python 2 — 2.7. Мы будем использовать Python 3. Отличия в синтаксисе между 2-й и 3-й версиями минимальны, но есть. С точки зрения обучения, неважно, что использовать — переключиться с одного на другой — дело получаса.
  • Anaconda — это дистрибутив, включающий в себя Python и всякие другие полезные для нас штуки, в том числе множество научных библиотек, которыми мы будем пользоваться.
  • Jupyter — это среда, с помощью которой можно работать с Python в интерактивном режиме: ввели команду — получили результат. Раньше она называлась IPython Notebook, но сейчас она умеет работать с другими языками и её переименовали в Jupyter.
  • Блокнотами (notebook) называются документы, получающиеся при такой интерактивной работе. Они состоят из ячеек с кодом, результатами выполнения кода и дополнительными комментариями. Этот документ также является блокнотом. Блокноты имеют расширение ipynb.

Картинка для привлечения внимания

Подробнее о картинках см. ниже.

In [1]:
import matplotlib.pyplot as plt
%matplotlib inline 
%config InlineBackend.figure_format = 'svg'
import numpy as np
r = np.linspace(0,200,100)
ax = plt.subplot(111, projection='polar')
ax.plot(r, r, color='b', linewidth=1.5);

Базовый Python

Этот раздел содержит шпаргалку по базовым возможностям Python. Вы можете его пропустить (сразу перейти к разделу numpy) и обращаться к нему по мере необходимости, или кратко просмотреть, чтобы узнать, чем Python отличается от других языков.

Арифметика

In [2]:
n = 10 # целое число (int)
x = 0.12 # число с плавающей точкой (float)
n * x # Результат выполнения последней команды будет выведем автоматически.
Out[2]:
1.2
In [3]:
type(n) # int — целое
Out[3]:
int
In [4]:
z = 1 + 2j # комплексное число, j — это мнимая единица
z ** 2 # возведение в степень — две звёздочки, а не крышечка; крышечка это побитный XOR
Out[4]:
(-3+4j)
In [5]:
print(5 / 2)  # деление в числах с плавающей точкой
print(5 // 2) # целочисленное деление
print(5 % 2)  # взятие остатка
2.5
2
1
In [6]:
from math import sqrt
# импорт одной функции из модуля (библиотеки)

sqrt(4)
Out[6]:
2.0
In [7]:
import math
# другой способ импорта

math.sin(math.pi)
Out[7]:
1.2246467991473532e-16
In [8]:
import math as m
# то же, что и предыдущий, но короче

m.sin(m.pi)
Out[8]:
1.2246467991473532e-16

Подробнее: конспект, видео.

Структуры данных

Списки

In [9]:
my_list = [4, 10, 1, 3]
# список (что-то вроде динамического массива)

print(my_list[2])
# нумерация начинается с нуля
1
In [10]:
my_list[3] = 99
my_list
Out[10]:
[4, 10, 1, 99]
In [11]:
my_list.append(100)
my_list
Out[11]:
[4, 10, 1, 99, 100]
In [12]:
my_list + [2, 3, 4] # конкатенация
Out[12]:
[4, 10, 1, 99, 100, 2, 3, 4]
In [13]:
my_list # не изменился
Out[13]:
[4, 10, 1, 99, 100]
In [14]:
my_list.extend([2, 3, 4])
my_list
Out[14]:
[4, 10, 1, 99, 100, 2, 3, 4]
In [15]:
my_list[2:5] # срез (slice); первый элемент включается, последний не включается
Out[15]:
[1, 99, 100]
In [16]:
min(my_list)
Out[16]:
1
In [17]:
max(my_list)
Out[17]:
100
In [18]:
sum(my_list)
Out[18]:
223
In [19]:
sorted(my_list)
Out[19]:
[1, 2, 3, 4, 4, 10, 99, 100]
In [20]:
my_list # список не изменился
Out[20]:
[4, 10, 1, 99, 100, 2, 3, 4]
In [21]:
my_list.sort() 
my_list # а теперь изменился
Out[21]:
[1, 2, 3, 4, 4, 10, 99, 100]

Подробнее: конспект, видео.

Кортежи

Кортеж — это неизменяемый список.

In [22]:
my_tuple = (12, 8, 3)
my_tuple[1]
Out[22]:
8
In [23]:
my_tuple[1] = 10
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-23-78e14473086f> in <module>()
----> 1 my_tuple[1] = 10

TypeError: 'tuple' object does not support item assignment

Словари

In [24]:
my_dict = {'Alice': 5, 'Bob': 3, 'Claudia': 4}
my_dict['Alice']
Out[24]:
5
In [25]:
my_dict['Bob'] = 4
my_dict
Out[25]:
{'Alice': 5, 'Bob': 4, 'Claudia': 4}

Подробнее: конспект, видео.

Множества

In [26]:
A = {9, 2, 5, 3, 10}
B = {2, 10, 12, 15}
9 in A
Out[26]:
True
In [27]:
A | B # объединение
Out[27]:
{2, 3, 5, 9, 10, 12, 15}
In [28]:
A & B # пересечение
Out[28]:
{2, 10}

Управляющие конструкции

Проверка условий

In [29]:
x = 10
if x > 8:
    print("x is rather big")
else:
    print("x is very small")
    print("Could you increase it?")
print("x =", x)
x is rather big
x = 10

Блок, относящийся к управляющей конструкции, выделяется отступом. Закончился отступ — закончился и блок. Никаких фигурных скобок или конструкций begin-end в Python не предусмотрено.

In [30]:
x = 15
if x > 10 and x % 2 == 0:
    print("x is big and even")
else:
    print("either x is not so big or it is not even")
either x is not so big or it is not even

Подробнее: конспект.

Циклы

Цикл for

In [31]:
my_list = [8, 9, 12]
for x in my_list:
    print(x)
    print("Let's go to the next x")
8
Let's go to the next x
9
Let's go to the next x
12
Let's go to the next x
In [32]:
for x in range(5):
    print(x)
0
1
2
3
4
In [33]:
list(range(2, 8))
Out[33]:
[2, 3, 4, 5, 6, 7]
In [34]:
# посчитаем факториал
s = 1
n = 15
for x in range(1, n + 1):
    s = s * x
print(s)
1307674368000
In [35]:
# хотя конечно проще так:
from math import factorial
factorial(15)
Out[35]:
1307674368000
In [36]:
for i, x in enumerate(my_list):
    print("my_list[{0}] = {1}".format(i, x))
my_list[0] = 8
my_list[1] = 9
my_list[2] = 12
In [37]:
other_list = [10, 12, 13]
# zip — застёжка-молния, состёгивает два или несколько списков
for x, y in zip(my_list, other_list):
    print("x = {}, y = {}".format(x, y))
    print("x + y =", x + y)
x = 8, y = 10
x + y = 18
x = 9, y = 12
x + y = 21
x = 12, y = 13
x + y = 25

Цикл while

In [38]:
# найдем последнее число Фибоначчи, меньшее 1000
a = 1
b = 1
while b < 1000:
    c = a + b
    a = b
    b = c
    # можно было бы так: a, b = b, a + b
print(a)
987
In [39]:
# найдем предпоследнее число Фибоначчи, меньшее 1000
a = 1
b = 1
while True: # выполнять всегда
    c = a + b
    if c > 1000:
        break
    a = b
    b = c
print(a)
610

Подробнее: while.

Математические библиотеки

numpy

numpy — это библиотека для эффективной работы с массивами. Массивы numpy называются numpy.array или numpy.ndarray (это почти одно и то же) и похожи на списки, но работают быстрее. По своей концепции они похожи на MATLAB.

In [40]:
import numpy as np
x = np.array([4, 3, 10, 3])
In [41]:
x[2]
Out[41]:
10
In [42]:
y = np.array([10, 12, 3, 15])
x + y # поэлементное сложение
Out[42]:
array([14, 15, 13, 18])
In [43]:
x * y # и умножение тоже
Out[43]:
array([40, 36, 30, 45])
In [44]:
x.dot(y) # скалярное произведение
Out[44]:
151
In [45]:
x @ y # так тоже можно
Out[45]:
151
In [46]:
A = np.array([[1, 2], [3, 4]]) # матрица
u = np.array([[1, 0]])
A @ u.T # умножение матрицы на вектор
Out[46]:
array([[1],
       [3]])
In [47]:
math.sin(x)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-47-36f3cd9cf9e1> in <module>()
----> 1 math.sin(x)

TypeError: only length-1 arrays can be converted to Python scalars
In [48]:
np.sin(x) # поэлементное применение
Out[48]:
array([-0.7568025 ,  0.14112001, -0.54402111,  0.14112001])
In [49]:
x
Out[49]:
array([ 4,  3, 10,  3])
In [50]:
x[x<10] # выбрать все элементы x, меньшие 10
Out[50]:
array([4, 3, 3])

numpy позволяет в большинстве случаев обходиться без циклов (которые в Python довольно медленные), заменяя их на поэлементные операции над списками.

Подробнее: конспект, видео.

scipy

Всякая разная математика.

Линейная алгебра

Решим систему уравнений

\begin{equation} \begin{pmatrix} 3 & 4\\ 4 & -3 \end{pmatrix} \begin{pmatrix} x\ y

\end{pmatrix}
=
\begin{pmatrix}
10\\
5
\end{pmatrix}

\end{equation}

In [51]:
from scipy.linalg import solve
solve(np.array([[3, 4], [4, -3]]), np.array([10, 5]).T) # решить систему линейных уравнений
Out[51]:
array([ 2.,  1.])

Подробнее: документация (англ.)

Интегрирование и решение дифференциальных уравнений

Посчитаем $\int_0^1 x^2 \; dx$.

In [52]:
from scipy.integrate import quad
def f(x): # так мы определили функцию f(x)=x**2
    return x**2
quad(f, 0, 1)
Out[52]:
(0.33333333333333337, 3.700743415417189e-15)

Решим дифференциальное уравнение $\dot x = x$.

In [53]:
from scipy.integrate import odeint
T = np.linspace(0, 1, 10) # равномерное разбиение отрезка [0, 1] с помощью 10 точек
def f(x, t):
    return x
odeint(f, 1, T)
Out[53]:
array([[ 1.        ],
       [ 1.11751906],
       [ 1.24884886],
       [ 1.39561243],
       [ 1.55962349],
       [ 1.742909  ],
       [ 1.94773405],
       [ 2.17663003],
       [ 2.43242551],
       [ 2.7182819 ]])

Подробнее: документация (англ.)

Решим систему \begin{equation} \dot x = y,\quad \dot y = - x \end{equation}

In [54]:
def f(X, t):
    # x = X[0], y = X[1]
    return [X[1], -X[0]]
T = np.linspace(0, 2 * m.pi, 10)
odeint(f, [1, 0], T)
Out[54]:
array([[  1.00000000e+00,   0.00000000e+00],
       [  7.66044451e-01,  -6.42787586e-01],
       [  1.73648228e-01,  -9.84807729e-01],
       [ -4.99999898e-01,  -8.66025431e-01],
       [ -9.39692552e-01,  -3.42020255e-01],
       [ -9.39692654e-01,   3.42020020e-01],
       [ -5.00000115e-01,   8.66025338e-01],
       [  1.73648045e-01,   9.84807790e-01],
       [  7.66044373e-01,   6.42787736e-01],
       [  1.00000004e+00,   1.47816225e-07]])

sympy

Символьная математика.

In [55]:
import sympy as sp
from sympy import init_printing
init_printing(use_latex = 'mathjax')
In [56]:
x, y, z = sp.symbols('x y z')
sp.expand((x + y)**10)
Out[56]:
$$x^{10} + 10 x^{9} y + 45 x^{8} y^{2} + 120 x^{7} y^{3} + 210 x^{6} y^{4} + 252 x^{5} y^{5} + 210 x^{4} y^{6} + 120 x^{3} y^{7} + 45 x^{2} y^{8} + 10 x y^{9} + y^{10}$$
In [57]:
sp.exp(-x**2 / 2).diff(x)
Out[57]:
$$- x e^{- \frac{x^{2}}{2}}$$
In [58]:
sp.sin(x).series()
Out[58]:
$$x - \frac{x^{3}}{6} + \frac{x^{5}}{120} + \mathcal{O}\left(x^{6}\right)$$
In [59]:
sp.exp(x**2 + x**3).series()
Out[59]:
$$1 + x^{2} + x^{3} + \frac{x^{4}}{2} + x^{5} + \mathcal{O}\left(x^{6}\right)$$
In [60]:
sp.integrate(sp.sin(x) ** 10)
Out[60]:
$$\frac{63 x}{256} - \frac{1}{10} \sin^{9}{\left (x \right )} \cos{\left (x \right )} - \frac{9}{80} \sin^{7}{\left (x \right )} \cos{\left (x \right )} - \frac{21}{160} \sin^{5}{\left (x \right )} \cos{\left (x \right )} - \frac{21}{128} \sin^{3}{\left (x \right )} \cos{\left (x \right )} - \frac{63}{256} \sin{\left (x \right )} \cos{\left (x \right )}$$
In [61]:
sp.exp(x**10).diff(x).subs({x:1})
Out[61]:
$$10 e$$
In [62]:
print(sp.exp(x**10).diff(x).subs({x:1}).evalf(n=1000))
27.18281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702761838606261331384583000752044933826560297606737113200709328709127443747047230696977209310141692836819025515108657463772111252389784425056953696770785449969967946864454905987931636889230098793127736178215424999229576351482208269895193668033182528869398496465105820939239829488793320362509443117301238197068416140397019837679320683282376464804295311802328782509819455815301756717361332069811250996181881593041690351598888519345807273866738589422879228499892086805825749279610484198444363463244968487560233624827041978623209002160990235304369941849146314093431738143640546253152096183690888707016768396424378140592714563549061303107208510383750510115747704171898610687396965521267154688957035035
In [63]:
M = sp.Matrix([[x, 1], [1, x]])
M
Out[63]:
$$\left[\begin{matrix}x & 1\\1 & x\end{matrix}\right]$$
In [64]:
sp.Matrix([[x, 1], [1, x]]).eigenvects()
Out[64]:
$$\left [ \left ( x - 1, \quad 1, \quad \left [ \left[\begin{matrix}-1\\1\end{matrix}\right]\right ]\right ), \quad \left ( x + 1, \quad 1, \quad \left [ \left[\begin{matrix}1\\1\end{matrix}\right]\right ]\right )\right ]$$

Подробнее: документация (англ.)

matplotlib

Картинки. Основная библиотека — matplotlib.pyplot — очень похожа на MATLAB и работает в тесной связке с numpy и другими библиотеками.

In [65]:
%matplotlib inline 
# иначе картинки не будут рисоваться
%config InlineBackend.figure_format = 'svg'
# вывод в SVG для пущей красивости
import matplotlib.pyplot as plt
plt.plot([0, 1, 2, 3], [0, 1, 4, 9])
Out[65]:
[<matplotlib.lines.Line2D at 0x108575b00>]
In [66]:
x = np.linspace(-4, 4)
plt.plot(x, x**2)
Out[66]:
[<matplotlib.lines.Line2D at 0x10867c4e0>]
In [67]:
plt.plot(x, np.sin(x**2))
Out[67]:
[<matplotlib.lines.Line2D at 0x10877be80>]
In [68]:
x = np.linspace(-4, 4, 1000) # возьмём точек побольше
plt.plot(x, np.sin(x**2))
Out[68]:
[<matplotlib.lines.Line2D at 0x10887f198>]
In [69]:
plt.plot(np.sin(x))
plt.plot(np.cos(x))
Out[69]:
[<matplotlib.lines.Line2D at 0x1088894a8>]
In [70]:
plt.plot(x, 1/x)
Out[70]:
[<matplotlib.lines.Line2D at 0x108a6bdd8>]
In [71]:
plt.ylim(-4, 4)
plt.plot(x, 1/x, label='$y = 1/x$')
plt.plot(x, np.zeros_like(x))
plt.legend()
Out[71]:
<matplotlib.legend.Legend at 0x108b70908>
In [72]:
plt.figure(figsize=(6, 6)) # квадратненько
x = np.linspace(-1.5, 1.5, 100)
y = x**2
plt.plot(x, y, label='$y = x^2$')
plt.plot(y, x, label='$x = y^2$')
plt.xlabel("$x$", fontsize=20)
plt.ylabel("$y$", fontsize=20)
plt.legend()
Out[72]:
<matplotlib.legend.Legend at 0x108d16908>

Scatter plot

In [73]:
from scipy.stats import norm
x = norm.rvs(size = 500)
y = x + norm.rvs(size = 500)
plt.scatter(x, y)
Out[73]:
<matplotlib.collections.PathCollection at 0x108c6b1d0>

Полярные координаты

In [74]:
r = np.linspace(0, 3.0, 10000)
theta = 2 * np.pi * np.exp(r)

ax = plt.subplot(111, projection='polar')
ax.plot(theta, r, color='r', linewidth=1)
Out[74]:
[<matplotlib.lines.Line2D at 0x10936b438>]

Векторные поля и ОДУ

Векторное поле дифференциального уравнения:

$$\dot x=-0.1x+y,\quad \dot y = -x - 0.1y$$

In [75]:
plt.figure(figsize=(6,6))
plt.xlabel("$x$", fontsize=20)
plt.ylabel("$y$", fontsize=20)
x, y = np.mgrid[-3:3:0.5, -3:3:0.5]
plt.quiver(x, y, -0.1 * x + y, -x - 0.1 * y)
Out[75]:
<matplotlib.quiver.Quiver at 0x109378860>

Фазовый портрет того же уравнения.

In [85]:
plt.figure(figsize=(6,6))
plt.xlabel("$x$", fontsize=20)
plt.ylabel("$y$", fontsize=20)

y, x = np.mgrid[-3:3:21j, -3:3:21j] # 21j здесь означает, что нужна 21 точка

# Обратите внимание: y и x идут в обратном порядке!

plt.streamplot(x, y, -0.1 * x + y, -x - 0.1 * y)
Out[85]:
<matplotlib.streamplot.StreamplotSet at 0x10f9e8a58>

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

In [76]:
def f(X, t):
    return np.array([[-0.1, 1], [-1, -0.1]] @ X)
T = np.linspace(0, 10*m.pi, 1000)
inits = [[1, 0], [0, 1], [1, 1]]
plt.figure(figsize=(6,6))
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xlabel("$x$", fontsize=20)
plt.ylabel("$y$", fontsize=20)
for init in inits:
    traj = odeint(f, init, T)
    plt.plot(traj[:, 0], traj[:, 1])

Линии уровня (contour plot)

Линии уровня функции $z=xy$.

In [77]:
x, y = np.mgrid[-3:3:0.01, -3:3:0.01]
z = x * y
plt.contour(x, y, z, 20, cmap='gnuplot')
Out[77]:
<matplotlib.contour.QuadContourSet at 0x10a42df28>

Касп (он же «ласточкин хвост») $x^2 - y^3 = 0$.

In [78]:
import matplotlib
matplotlib.rc('font', family='Arial') # иначе русские буквы не сработают
x, y = np.mgrid[-1:1:0.001, 0:1:0.001]
z = x**2 - y**3
plt.contour(x, y, z, levels=[0])
plt.title("Ласточкин хвост")
Out[78]:
<matplotlib.text.Text at 0x10a73e710>

Подробнее: pyplot tutorial (англ.)

Виджеты

Виджеты позволяют добавить интерактивности: создают элементы управления (например, слайдеры), с помощью которых вы можете настраивать параметры вызываемых функций и сразу получать результат. Виджеты работают только если вы загрузили блокнот и свой IPython Notebook или в tmpnb. Если вы просто читаете эту страничку в Интернете, скорее всего, вы ничего интересного не увидите.

In [1]:
import matplotlib.pyplot as plt
import numpy as np
In [2]:
%matplotlib inline
In [3]:
from ipywidgets import interact, interactive, fixed, FloatSlider
import ipywidgets as widgets
In [4]:
@interact(a=FloatSlider(min=0, max=10, value=1, step=1e-3), 
          b=FloatSlider(min=0, max=10, value=1, step=1e-3))
def plot_sin(a, b):
    x = np.linspace(-4,4,300)
    plt.ylim(-2,2)
    plt.plot(np.sin(a*x)*b)

Подробнее: виджеты (англ.)