Symbolický a sympatický SymPy

Už jsme viděli mnoho nástrojů pro numerické výpočty. Co ale symbolické výpočty? Každý ví, že derivovat umí i cvičená opice, bude to umět i Python? Nečekaná odpověď je ano. Symbolické výpočty naučil Python Ondřej Čertík v balíku SymPy.

Tento notebook byl z (velké) části převzat a přeložen z J.R. Johansson: Lectures on scientific computing with Python - díky.

Na úvod

Někteří z vás možná znáte nějaký systém pro počítačovou algebru (Computer Algebra Systems -- CAS), např. Maple, Mathematica, Derive, Maxima, Reduce. Pro Python existují dva velké projekty počítačové algebry:

  • SymPy - modul který může být použit v jakémkoli Python programu a je dobře podporován v IPython Notebook.
  • Sage - toto je už kompletní (a velice obsáhlý) systém, který si klade za cíl být open source konkurentem komerčním produktům.

My se tady podíváme na některé základní možnosti SymPy.

In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
In [2]:
from sympy import *
In [3]:
# zastaralé
# %load_ext sympy.interactive.ipythonprinting
from sympy import init_printing
init_printing()

Definujeme symboly

Pro symbolické výpočty potřebujeme pochopitelně symboly, tak jak jsme zvyklí už z matematiky na základní škole. V Pythonu samotném máme sice proměnné, které jsou v podstatě také symboly, ale operace s nimy se řídí zcela jinými pravidly než potřebujeme pro symbolické výpočty. Naštěstí tu je třída sympy.Symbol.

In [4]:
x = Symbol('x')
x
Out[4]:
$$x$$

Co když napíšeme něco trochu složitějšího.

In [5]:
(pi + x) / 2
Out[5]:
$$\frac{x}{2} + \frac{\pi}{2}$$
In [6]:
# co jsme to vůbec dostali za typ
type(_)
Out[6]:
sympy.core.add.Add

Můžeme také přičadit symbolům nějaké vlastnosti (to se pak pochopitelně může projevit v dalších výpočtech).

In [7]:
a = Symbol('a', real=True)
a.is_complex
Out[7]:
True
In [8]:
b = Symbol('b', positive=True)
In [9]:
b.is_real
Out[9]:
True
In [10]:
b > 0
Out[10]:
$$\mathrm{True}$$

Zlomky

In [11]:
r1 = Rational(4,5)
r2 = Rational(5,4)
r1, r2
Out[11]:
$$\left ( \frac{4}{5}, \quad \frac{5}{4}\right )$$
In [12]:
r1 + r2
Out[12]:
$$\frac{41}{20}$$

Vyčíslování

In [13]:
y = (x + pi)**2

Numerickou hodnotu můžeme získat pomocí funkce N. Často také využijeme metodu subs:

In [14]:
N(y.subs(x, 2), 5)
Out[14]:
$$26.436$$

To samé pomocí metody evalf. Pro obojí můžeme zadat počet platných číslic.

In [15]:
pi.evalf(100)
Out[15]:
$$3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068$$

Pokud chceme vytvořit ze symbolického výrazu funkci, použijeme lambdify:

In [16]:
# první argument je seznam proměnných (podobně jako pro lambda funkce)
f_x = lambdify([x], (x + pi)**2, 'numpy')
In [17]:
xa = linspace(-10, 10)
fix, ax = subplots()
ax.plot(xa, f_x(xa))
Out[17]:
[<matplotlib.lines.Line2D at 0x7f640e3ca6d8>]

Symbolické úpravy

Toto je velice důležitá aplikace, která nám může v mnoha případech ušetřit nemálo práce.

Expand a factor

Začněme pracovat s polynomem, zadaným jako

In [18]:
y = (x+1)*(x+2)*(x+3)
y
Out[18]:
$$\left(x + 1\right) \left(x + 2\right) \left(x + 3\right)$$

Polynom rozvineme pomocí expand:

In [19]:
z = expand(y)
z
Out[19]:
$$x^{3} + 6 x^{2} + 11 x + 6$$

Pomocí factor můžeme dostat zpět původní faktorizovaný výraz.

In [20]:
factor(z)
Out[20]:
$$\left(x + 1\right) \left(x + 2\right) \left(x + 3\right)$$

expand můžeme použít i pro trigonometrické funkce:

In [21]:
expand(sin(a+b), trig=True)
Out[21]:
$$\sin{\left (a \right )} \cos{\left (b \right )} + \sin{\left (b \right )} \cos{\left (a \right )}$$

Zjednodušování pomocí simplify

In [22]:
# tohle by měla být hračka
simplify(sin(a)**2 + cos(a)**2)
Out[22]:
$$1$$

Derivace a integrály

SymPy umí symbolicky derivovat (je tedy aspoň tak dobrý jako cvičená opice) a i integrovat.

In [23]:
y = (x**2 + sin(x))**2
y
Out[23]:
$$\left(x^{2} + \sin{\left (x \right )}\right)^{2}$$
In [24]:
diff(y, x)
Out[24]:
$$\left(4 x + 2 \cos{\left (x \right )}\right) \left(x^{2} + \sin{\left (x \right )}\right)$$

Derivovat můžeme i funkce více proměnných.

In [25]:
x = Symbol('x')
y = Symbol('y')
z = cos(y) * (x**3 + 2*x**2*y)
z
Out[25]:
$$\left(x^{3} + 2 x^{2} y\right) \cos{\left (y \right )}$$

Tohle spočítá

$\displaystyle \frac{{{{\rm{d}}^3}z}}{{{\rm{d}}x{\rm{d}}{y^2}}} $

In [26]:
diff(z, x, 1, y, 2)
Out[26]:
$$- x \left(\left(3 x + 4 y\right) \cos{\left (y \right )} + 8 \sin{\left (y \right )}\right)$$

Integrace

In [27]:
f = sin(x*y)*cos(x)
f
Out[27]:
$$\sin{\left (x y \right )} \cos{\left (x \right )}$$
In [30]:
integrate(f, x)
Out[30]:
$$\begin{cases} \frac{1}{2} \cos^{2}{\left (x \right )} & \text{for}\: y = -1 \\- \frac{1}{2} \cos^{2}{\left (x \right )} & \text{for}\: y = 1 \\- \frac{y \cos{\left (x \right )}}{y^{2} - 1} \cos{\left (x y \right )} - \frac{\sin{\left (x \right )}}{y^{2} - 1} \sin{\left (x y \right )} & \text{otherwise} \end{cases}$$
In [31]:
integrate(f, y)
Out[31]:
$$\left(\begin{cases} 0 & \text{for}\: x = 0 \\- \frac{1}{x} \cos{\left (x y \right )} & \text{otherwise} \end{cases}\right) \cos{\left (x \right )}$$

Generace kódu

Automatická generace kódu vlastnost, kterou oceníme ve chvíli, kdy cheme implementovat naše analytické výsledky v numerické simulaci. Místo abychom začali ručně přepisovat do programovacího jazyka jako je např. Fortran nebo C, může SymPy tuto nezábavnou práci udělat za nás. Navíc při tom neudělá chyby (pravděpodobně).

Pozor, tento module je work in progress.

In [32]:
# řekněme že chceme někde použít tento výsledek
f = sin(x*y**2)*exp(y)
f
Out[32]:
$$e^{y} \sin{\left (x y^{2} \right )}$$
In [33]:
from sympy.utilities.codegen import codegen
In [34]:
f_source = codegen(("f_fortran", f), "F95", "f_fortran")
print(f_source[0][1])
!******************************************************************************
!*                      Code generated with sympy 0.7.6                       *
!*                                                                            *
!*              See http://www.sympy.org/ for more information.               *
!*                                                                            *
!*                       This file is part of 'project'                       *
!******************************************************************************

REAL*8 function f_fortran(x, y)
implicit none
REAL*8, intent(in) :: x
REAL*8, intent(in) :: y

f_fortran = exp(y)*sin(x*y**2)

end function

In [35]:
f_source = codegen(("f_fortran", f), "C", "f_fortran")
print(f_source[0][1])
/******************************************************************************
 *                      Code generated with sympy 0.7.6                       *
 *                                                                            *
 *              See http://www.sympy.org/ for more information.               *
 *                                                                            *
 *                       This file is part of 'project'                       *
 ******************************************************************************/
#include "f_fortran.h"
#include <math.h>

double f_fortran(double x, double y) {

   double f_fortran_result;
   f_fortran_result = exp(y)*sin(x*pow(y, 2));
   return f_fortran_result;

}

Další možnosti SymPy

Ukázali jsme si základy práce se symbolickými výpočty pomocí SymPy. Není v našich silách ukázat, co všechno SymPy umí -- je toho opravdu hodně. Přehled můžeme získat, podívame-li se na obsah v dokumentaci:

  • SymPy Core
  • Combinatorics Module
  • Number Theory
  • Concrete Mathematics
  • Numerical evaluation
  • Functions Module
  • Geometry Module
  • Geometric Algebra Module
  • Geometric Algebra Module for SymPy
  • Extended LaTeXModule for SymPy
  • Symbolic Integrals
  • Numeric Integrals
  • Logic Module
  • Matrices
  • Mpmath
  • Polynomials Manipulation Module
  • Printing System
  • Plotting Module
  • Pyglet Plotting Module
  • Assumptions module
  • Term rewriting
  • Series Expansions
  • Sets
  • Simplify
  • Details on the Hypergeometric Function Expansion Module
  • Statistics
  • Stats
  • ODE
  • PDE
  • Solvers
  • Tensor Module
  • Utilities
  • Parsing input
  • Physics Module
  • Category Theory Module
  • Differential Geometry Module
  • Contributions to docs

Komentáře

Comments powered by Disqus