Программирование и научные вычисления на языке Python/§19

Хотя численные методы, некоторые из которых описаны в предыдущих уроках, и находят широкое применение, большей точностью и наглядностью обладают символьные вычисления, которые работают с математическими выражениями, собственно, как и предполагает математика, как с последовательностями символов. Системы, занимающиеся символьными вычислениями, называют также системами компьютерной алгебры. Примерами таких систем служат известные математические среды Maple, Mathcad, Mathematica, Maxima и т. д. В качестве инструмента символьных вычислений мы рассмотрим библиотеку SymPy.

SymPy представляет собой библиотеку символьных вычислений, которая в конечной цели стремится стать полноценной системой компьютерной алгебры, сохраняя при этом как можно более простой код, ясный для понимания и дальнейших изменений и дополнений. SymPy написан исключительно на Python и не требует никаких других библиотек. Установка SymPy происходит аналогично другим продуктам, рассмотренным ранее.

isympy

Если вы уже используете для своих научных экспериментов IPython, то для вас окажется приятным тот факт, что в SymPy для него имеется обертка — isympy, запуск которой в консоли позволяет упростить некоторые стандартные действия относительно переменных и представлять данные на выходе в удобном для понимания формате записи:

leo@leo:~$ isympy
Python 2.6.5 console for SymPy 0.6.6

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z = symbols('xyz')
>>> k, m, n = symbols('kmn', integer=True)
>>> f, g, h = map(Function, 'fgh')

Documentation can be found at http://sympy.org/

In [1]: (1/cos(x)).series(x, 0, 10)
Out[1]: 
     2      4       6        8           
    x    5⋅x    61⋅x    277⋅x            
1 + ── + ──── + ───── + ────── + O(x**10)
    2     24     720     8064

Далее в уроке, в зависимости от удобства представления, одни данные будут показаны как при работе в стандартном IDLE, другие — как при работе в isympy.

Калькулятор

SymPy содержит три встроенных числовых типа: Real, Rational и Integer. Класс Rational отвечает за дроби, любую пару целых чисел он определяет как числитель и знаменатель, далее над дробями можно производить типичные действия: умножать, делить, складывать и так далее, при этом например:

>>> from sympy import *
>>> a = Rational(1,2)

>>> a
1/2

>>> a*2
1

>>> Rational(2)**50/Rational(10)**50
1/88817841970012523233890533447265625


# при этом:
>>> 1/2
0
>>> 1.0/2
0.5


# можно использовать математические константы
>>> pi**2
pi**2

>>> pi.evalf()
3.14159265358979

>>> (pi+exp(1)).evalf()
5.85987448204884


# а также и символ для бесконечности
>>> oo > 99999
True
>>> oo + 1
oo

Символы

В отличие от других систем компьютерной алгебры, в SymPy необходимо предварительно обозначить, какие переменные являются символьными:

>>> from sympy import *
>>> x = Symbol('x')
>>> y = Symbol('y')

Теперь вы можете играть с ними:

>>> x+y+x-y
2*x

>>> (x+y)**2
(x + y)**2

>>> ((x+y)**2).expand()
2*x*y + x**2 + y**2

И заменять их другими символами или числами, используя конструкцию subs(old, new):

>>> ((x+y)**2).subs(x, 1)
(1 + y)**2

>>> ((x+y)**2).subs(x, y)
4*y**2

Разложение на элементарные

Для разложения дроби на элементарные используется apart(expr, x):

In [1]: 1/( (x+2)*(x+1) )
Out[1]:
       1
───────────────
(2 + x)*(1 + x)

In [2]: apart(1/( (x+2)*(x+1) ), x)
Out[2]:
  1       1
───── - ─────
1 + x   2 + x

In [3]: (x+1)/(x-1)
Out[3]:
-(1 + x)
────────
 1 - x

In [4]: apart((x+1)/(x-1), x)
Out[4]:
      2
1 - ─────
    1 - x

Чтобы снова представить их в виде одной дроби together(expr, x):

In [7]: together(1/x + 1/y + 1/z)
Out[7]:
x*y + x*z + y*z
───────────────
     x*y*z

In [8]: together(apart((x+1)/(x-1), x), x)
Out[8]:
-1 - x
──────
1 - x

In [9]: together(apart(1/( (x+2)*(x+1) ), x), x)
Out[9]:
       1
───────────────
(2 + x)*(1 + x)

Пределы

Пределы в SymPy определяются очень просто. Чтобы найти предел функции f(x) при x -> 0, необходимо так и записать limit(f(x), x, 0):

>>> from sympy import *
>>> x=Symbol("x")
>>> limit(sin(x)/x, x, 0)
1

Естественно, можно использовать и предел в бесконечности:

>>> limit(x, x, oo)
oo

>>> limit(1/x, x, oo)
0

>>> limit(x**x, x, 0)
1

Дифференцирование

Любое SymPy выражение можно продифференцировать, используя diff(func, var):

>>> from sympy import *
>>> x = Symbol('x')
>>> diff(sin(x), x)
cos(x)
>>> diff(sin(2*x), x)
2*cos(2*x)

>>> diff(tan(x), x)
1 + tan(x)**2

По определению производной вы можете проверить правильность решения:

>>> limit((tan(x+y)-tan(x))/y, y, 0)
1 + tan(x)**2

Вторые, третьи и так далее производные можно найти с помощью diff(func, var, n):

>>> diff(sin(2*x), x, 1)
2*cos(2*x)

>>> diff(sin(2*x), x, 2)
-4*sin(2*x)

>>> diff(sin(2*x), x, 3)
-8*cos(2*x)

Разложение в ряд

Используем метод .series(var, point, order):

>>> from sympy import *
>>> x = Symbol('x')
>>> cos(x).series(x, 0, 10)
1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320 + O(x**10)
>>> (1/cos(x)).series(x, 0, 10)
1 + x**2/2 + 5*x**4/24 + 61*x**6/720 + 277*x**8/8064 + O(x**10)

Другой простой пример:

from sympy import Symbol, pprint

x = Symbol("x")
y = Symbol("y")

e = 1/(x + y)
s = e.series(x, 0, 5)

print(s)
pprint(s)

дает нам после запуска:

1/y + x**2*y**(-3) + x**4*y**(-5) - x*y**(-2) - x**3*y**(-4) + O(x**5)
     2    4         3
1   x    x    x    x
─ + ── + ── - ── - ── + O(x**5)
y    3    5    2    4
    y    y    y    y

Интегрирование

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

>>> from sympy import *
>>> x, y = symbols('xy')

Интегралы от элементарных функций:

>>> integrate(6*x**5, x)
x**6
>>> integrate(sin(x), x)
-cos(x)
>>> integrate(log(x), x)
-x + x*log(x)
>>> integrate(2*x + sinh(x), x)
cosh(x) + x**2

И от специальных:

>>> integrate(exp(-x**2)*erf(x), x)
pi**(1/2)*erf(x)**2/4

Можно найти и определенный интеграл:

>>> integrate(x**3, (x, -1, 1))
0
>>> integrate(sin(x), (x, 0, pi/2))
1
>>> integrate(cos(x), (x, -pi/2, pi/2))
2

Или несобственный интеграл:

>>> integrate(exp(-x), (x, 0, oo))
1
>>> integrate(log(x), (x, 0, 1))
-1

Комплексные числа

>>> from sympy import Symbol, exp, I
>>> x = Symbol("x")
>>> exp(I*x).expand()
exp(I*x)
>>> exp(I*x).expand(complex=True)
I*exp(-im(x))*sin(re(x)) + cos(re(x))*exp(-im(x))
>>> x = Symbol("x", real=True)
>>> exp(I*x).expand(complex=True)
I*sin(x) + cos(x)

Функции

Тригонометрические

In [1]: sin(x+y).expand(trig=True)
Out[1]: cos(x)*sin(y) + cos(y)*sin(x)

In [2]: cos(x+y).expand(trig=True)
Out[2]: cos(x)*cos(y) - sin(x)*sin(y)

In [3]: sin(I*x)
Out[3]: ⅈ*sinh(x)

In [4]: sinh(I*x)
Out[4]: ⅈ*sin(x)

In [5]: asinh(I)
Out[5]:
π*ⅈ
───
 2

In [6]: asinh(I*x)
Out[6]: ⅈ*asin(x)

In [15]: sin(x).series(x, 0, 10)
Out[15]:
     3     5     7       9
    x     x     x       x
x - ── + ─── - ──── + ────── + O(x**10)
    6    120   5040   362880

In [16]: sinh(x).series(x, 0, 10)
Out[16]:
     3     5     7       9
    x     x     x       x
x + ── + ─── + ──── + ────── + O(x**10)
    6    120   5040   362880

In [17]: asin(x).series(x, 0, 10)
Out[17]:
     3      5      7       9
    x    3*x    5*x    35*x
x + ── + ──── + ──── + ───── + O(x**10)
    6     40    112     1152

In [18]: asinh(x).series(x, 0, 10)
Out[18]:
     3      5      7       9
    x    3*x    5*x    35*x
x - ── + ──── - ──── + ───── + O(x**10)
    6     40    112     1152

Сферические гармоники

In [1]: from sympy.abc import theta, phi

In [2]: Ylm(1, 0, theta, phi)
Out[2]:
  ⎽⎽⎽
╲╱ 3 *cos(θ)
────────────
      ⎽⎽⎽
  2*╲╱ π

In [3]: Ylm(1, 1, theta, phi)
Out[3]:
   ⎽⎽⎽           ⅈ*φ
-╲╱ 6 *│sin(θ)│*ℯ
────────────────────
          ⎽⎽⎽
      4*╲╱ π

In [4]: Ylm(2, 1, theta, phi)
Out[4]:
   ⎽⎽⎽⎽                  ⅈ*φ
-╲╱ 30 *│sin(θ)│*cos(θ)*ℯ
────────────────────────────
              ⎽⎽⎽
          4*╲╱ π

Факториалы и гамма-функции

In [1]: x = Symbol("x")

In [2]: y = Symbol("y", integer=True)

In [3]: factorial(x)
Out[3]: Γ(1 + x)

In [4]: factorial(y)
Out[4]: y!

In [5]: factorial(x).series(x, 0, 3)
Out[5]:
                    2           2    2  2
                   x *EulerGamma    π *x
1 - x*EulerGamma + ────────────── + ───── + O(x**3)
                         2            12

Дзета-функция Римана

In [18]: zeta(4, x)
Out[18]: ζ(4, x)

In [19]: zeta(4, 1)
Out[19]:
 4
π
──
90

In [20]: zeta(4, 2)
Out[20]:
      4
     π
-1 + ──
     90

In [21]: zeta(4, 3)
Out[21]:
        4
  17   π
- ── + ──
  16   90

Многочлены

Полиномы Чебышева, Лежандра, Эрмита:

In [1]: chebyshevt(2, x)
Out[1]:
        2
-1 + 2*x

In [2]: chebyshevt(4, x)
Out[2]:
       2      4
1 - 8*x  + 8*x

In [3]: legendre(2, x)
Out[3]:
          2
       3*x
-1/2 + ────
        2

In [4]: legendre(8, x)
Out[4]:
           2         4         6         8
 35   315*x    3465*x    3003*x    6435*x
─── - ────── + ─────── - ─────── + ───────
128     32        64        32       128

In [5]: assoc_legendre(2, 1, x)
Out[5]:
        ⎽⎽⎽⎽⎽⎽⎽⎽
       ╱      2
-3*x*╲╱  1 - x

In [6]: assoc_legendre(2, 2, x)
Out[6]:
       2
3 - 3*x

In [7]: hermite(3, x)
Out[7]:
           3
-12*x + 8*x

Дифференциальные уравнения

In [1]: f(x).diff(x, x) + f(x)
Out[1]:
   2
  d
─────(f(x)) + f(x)
dx dx

In [2]: dsolve(f(x).diff(x, x) + f(x), f(x))
Out[2]: C₁*sin(x) + C₂*cos(x)

Алгебраические уравнения

In [3]: solve(x**4 - 1, x)
Out[3]: [ⅈ, 1, -1, -ⅈ]

In [4]: solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
Out[4]: {y: 1, x: -3}

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

Матрицы создаются как экземпляры специального класса:

>>> from sympy import Matrix
>>> Matrix([[1,0], [0,1]])
[1, 0]
[0, 1]

В них можно использовать не только числа, но и символы:

>>> x = Symbol('x')
>>> y = Symbol('y')
>>> A = Matrix([[1,x], [y,1]])
>>> A
[1, x]
[y, 1]

>>> A**2
[1 + x*y,     2*x]
[    2*y, 1 + x*y]

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

Сопоставление с образцом

Сопоставление с образцом (англ. Pattern matching) — метод анализа списков или других структур данных на наличие в них заданных образцов. для проведения сопоставления используется метод .match(), с использованием объекта класса Wild. Метод возвращает словарь с соответствующим решением:

>>> from sympy import *
>>> x = Symbol('x')
>>> p = Wild('p')
>>> (5*x**2).match(p*x**2)
{p_: 5}

>>> q = Wild('q')
>>> (x**2).match(p*x**q)
{p_: 1, q_: 2}

Если сопоставление не дает результатов, возвращается объект None:

>>> x = Symbol('x')
>>> p = Wild('p', exclude=[1,x])
>>> print (x+1).match(x+p) # 1 is excluded
None
>>> print (x+1).match(p+1) # x is excluded
None
>>> print (x+1).match(x+2+p) # -1 is not excluded
{p_: -1}

Печать

Существует несколько способов того как выражения могут быть выведены на экран.

Стандартный print

Привычный вывод str(expression):

>>> from sympy import Integral
>>> from sympy.abc import x
>>> print x**2
x**2
>>> print 1/x
1/x
>>> print Integral(x**2, x)
Integral(x**2, x)

Pretty printing

Формат, по которому работает isympy по умолчанию в виде псевдографического ascii-форматирования. Немного больше примеров здесь.

>>> from sympy import Integral, pprint
>>> from sympy.abc import x
>>> pprint(x**2)
 2
x
>>> pprint(1/x)
1
-
x
>>> pprint(Integral(x**2, x))
  /
 |
 |  2
 | x  dx
 |
/

Для того, чтобы сделать pprint по умолчанию в стандартном интерпретаторе, производим следующую процедуру:

>>> from sympy import *
>>> import sys
>>> sys.displayhook = pprint
>>> var("x")
x
>>> x**3/3
 3
x
--
3
>>> Integral(x**2, x)
  /
 |
 |  2
 | x  dx
 |
/

Python printing

>>> from sympy.printing.python import python
>>> from sympy import Integral
>>> from sympy.abc import x
>>> print python(x**2)
x = Symbol('x')
e = x**2
>>> print python(1/x)
x = Symbol('x')
e = 1/x
>>> print python(Integral(x**2, x))
x = Symbol('x')
e = Integral(x**2, x)

LaTeX printing

LaTeX — наиболее популярный набор макрорасширений системы компьютерной вёрстки TeX, который облегчает набор сложных документов.

>>> from sympy import Integral, latex
>>> from sympy.abc import x
>>> latex(x**2)
x^{2}
>>> latex(x**2, mode='inline')
$x^{2}$
>>> latex(x**2, mode='equation')
\begin{equation}x^{2}\end{equation}
>>> latex(x**2, mode='equation*')
\begin{equation*}x^{2}\end{equation*}
>>> latex(1/x)
\frac{1}{x}
>>> latex(Integral(x**2, x))
\int x^{2}\,dx

MathML

MathML (от англ. Mathematical Markup Language, язык математической разметки) — это приложение XML, используемое для представления математических символов и формул в документах WWW. MathML рекомендован математической группой W3C.

>>> from sympy.printing.mathml import mathml
>>> from sympy import Integral, latex
>>> from sympy.abc import x
>>> print mathml(x**2)
<apply><power/><ci>x</ci><cn>2</cn></apply>
>>> print mathml(1/x)
<apply><power/><ci>x</ci><cn>-1</cn></apply>

sympy.printing

Кроме того, как видно из рассмотрения python printing доступен модуль sympy.printing. Интересные функции, доступные в этом модуле:

• pretty(expr), pretty_print(expr), pprint(expr): возвращают expr в pretty-формате, как показано выше;
• latex(expr), print_latex(expr): возвращают expr в LaTeX представлении;
• mathml(expr), print_mathml(expr): то же для MathML;
• print_gtk(expr): выводит expr в Gtkmathview — GTK виджет, показывающий MathML код.

Что дальше

Этот урок является частичным переводом пособия для новичков. Если вам интересно, что еще может SymPy, то вы можете приступить к более полному руководству пользователя или пробежаться по модулям.

Ссылки

• Официальный сайт SymPy — отсюда можно скачать пакет для установки
• Tutorial — основа этого урока
• SymPy User’s Guide
• SymPy Modules Reference
• Собственное wiki