SymPy

SymPy is a package for symbolic calculations in python, similar to Mathematica. It works with expressions containing symbols.

from sympy import *
init_printing()

Symbols are basic bricks used to construct expressions. Each symbol has a name used for printing expressions. Objects of the class Symbol should be created and assigned to python variables in order to be used in expressions. The symbol name and the name of the variable to which this symbol is assigned are two independent things, and one may write abc=Symbol('xyz'). But then one has to write abc in input expressions, while SymPy will write xyz in output ones, producing unnecessary confusion. The python variable name should better be the same as the symbol name.

In languages specifically designed for symbolic calculations, such as Mathematica, if a variable to which nothing has been assigned is used, it automatically means a symbol with the same name. Python has not been designed for symbolic calculations. If you use a variable to which nothing has been assigned, you will get an error message. Symbol objects have to be created explicitly.

x=Symbol('x')

a=x**2-1
a

$$x^{2} - 1$$

Several symbols can be defined at once. The string is split at spaces.

y,z=symbols('y z')

Let's substitute $y+1$ for $x$.

a.subs(x,y+1)

$$\left(y + 1\right)^{2} - 1$$

Polynomials and rational functions

SymPy does not expand brackets automatically. The function expand is used for this.

a=(x+y-z)**6
a

$$\left(x + y - z\right)^{6}$$

a=expand(a)
a

$$x^{6} + 6 x^{5} y - 6 x^{5} z + 15 x^{4} y^{2} - 30 x^{4} y z + 15 x^{4} z^{2} + 20 x^{3} y^{3} - 60 x^{3} y^{2} z + 60 x^{3} y z^{2} - 20 x^{3} z^{3} + 15 x^{2} y^{4} - 60 x^{2} y^{3} z + 90 x^{2} y^{2} z^{2} - 60 x^{2} y z^{3} + 15 x^{2} z^{4} + 6 x y^{5} - 30 x y^{4} z + 60 x y^{3} z^{2} - 60 x y^{2} z^{3} + 30 x y z^{4} - 6 x z^{5} + y^{6} - 6 y^{5} z + 15 y^{4} z^{2} - 20 y^{3} z^{3} + 15 y^{2} z^{4} - 6 y z^{5} + z^{6}$$

Degree of the polynomial $a$ in $x$.

degree(a,x)

$$6$$

Let's collect terms with same power of $x$ together.

collect(a,x)

$$x^{6} + x^{5} \left(6 y - 6 z\right) + x^{4} \left(15 y^{2} - 30 y z + 15 z^{2}\right) + x^{3} \left(20 y^{3} - 60 y^{2} z + 60 y z^{2} - 20 z^{3}\right) + x^{2} \left(15 y^{4} - 60 y^{3} z + 90 y^{2} z^{2} - 60 y z^{3} + 15 z^{4}\right) + x \left(6 y^{5} - 30 y^{4} z + 60 y^{3} z^{2} - 60 y^{2} z^{3} + 30 y z^{4} - 6 z^{5}\right) + y^{6} - 6 y^{5} z + 15 y^{4} z^{2} - 20 y^{3} z^{3} + 15 y^{2} z^{4} - 6 y z^{5} + z^{6}$$

Any polynomial with integer coefficients can be factorized into polynomials with integer coefficients (which cannot be factorized further). There exist efficient algorithms to do this.

a=factor(a)
a

$$\left(x + y - z\right)^{6}$$

SymPy does not automatically cancel ratios of polynomials by their greatest common divisor. The function cancel is used for this.

a=(x**3-y**3)/(x**2-y**2)
a

$$\frac{x^{3} - y^{3}}{x^{2} - y^{2}}$$

cancel(a)

$$\frac{x^{2} + x y + y^{2}}{x + y}$$

SymPy does not automatically bring sums of rational expressions to common denominator. The function together is used for this.

a=y/(x-y)+x/(x+y)
a

$$\frac{x}{x + y} + \frac{y}{x - y}$$

together(a)

$$\frac{x \left(x - y\right) + y \left(x + y\right)}{\left(x - y\right) \left(x + y\right)}$$

The function simplify tries to rewrite an expression in a simplest form. This concept is not well defined (different forms may be considered simplest in different contexts), and there exists no algorithm for such simplification. The function simplify works heuristically, and it is not possible to guess in advance what simplifications it will try. It is very convenient in interactive sessions in order to check if it will succeed in rewriting an expression in some reasonably good form. But it is not desirable to use it in programs. There one should better use more specialized functions which perform well defined expression transformations.

simplify(a)

$$\frac{x^{2} + y^{2}}{x^{2} - y^{2}}$$

Partial fraction decomposition with respect to $x$.

apart(a,x)

$$- \frac{y}{x + y} + \frac{y}{x - y} + 1$$

Let's substitute some values for the symbils $x$ and $y$.

a=a.subs({x:1,y:2})
a

$$- \frac{5}{3}$$

And how much is it numerically?

a.n()

$$-1.66666666666667$$

Elementary functions

SymPy automatically applies simplifications of elementary functions which are correct everywhere.

sin(-x)

$$- \sin{\left (x \right )}$$

cos(pi/4),tan(5*pi/6)

$$\left ( \frac{\sqrt{2}}{2}, \quad - \frac{\sqrt{3}}{3}\right )$$

SymPy can work with floating point numbers having arbitrarily high precision. Here is $\pi$ with 100 significant digits.

pi.n(100)

$$3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068$$

E is the base of natural logarithms.

log(1),log(E)

$$\left ( 0, \quad 1\right )$$

exp(log(x)),log(exp(x))

$$\left ( x, \quad \log{\left (e^{x} \right )}\right )$$

Why not $x$? Try $x=2\pi i$.

sqrt(0)

$$0$$

sqrt(x)**4,sqrt(x**4)

$$\left ( x^{2}, \quad \sqrt{x^{4}}\right )$$

Why not $x^2$? Try $x=i$.

Symbols can have certain properties. E.g., they can be positive. Then SymPy can simplify square roots better.

p,q=symbols('p q',positive=True)
sqrt(p**2)

$$p$$

sqrt(12*x**2*y),sqrt(12*p**2*y)

$$\left ( 2 \sqrt{3} \sqrt{x^{2} y}, \quad 2 \sqrt{3} p \sqrt{y}\right )$$

Let the symbol $n$ be integer (I is the imaginary unit).

n=Symbol('n',integer=True)
exp(2*pi*I*n)

$$1$$

The method rewrite tries to rewrite an expression in terms of a given function.

cos(x).rewrite(exp),exp(I*x).rewrite(cos)

$$\left ( \frac{e^{i x}}{2} + \frac{1}{2} e^{- i x}, \quad i \sin{\left (x \right )} + \cos{\left (x \right )}\right )$$

asin(x).rewrite(log)

$$- i \log{\left (i x + \sqrt{- x^{2} + 1} \right )}$$

The function trigsimp tries to rewrite a trigonometric expression in a simplest form. In programs it is better to use more specialized functions.

trigsimp(2*sin(x)**2+3*cos(x)**2)

$$\cos^{2}{\left (x \right )} + 2$$

The function expand_trig expands sines and cosines of sums and multiple angles.

expand_trig(sin(x-y)),expand_trig(sin(2*x))

$$\left ( \sin{\left (x \right )} \cos{\left (y \right )} - \sin{\left (y \right )} \cos{\left (x \right )}, \quad 2 \sin{\left (x \right )} \cos{\left (x \right )}\right )$$

The inverse transformation, rewriting products and powers of sines and cosines into expressions linear in these functions, is needed more often. Suppose we work with a truncated Fourier series.

a1,a2,b1,b2=symbols('a1 a2 b1 b2')
a=a1*cos(x)+a2*cos(2*x)+b1*sin(x)+b2*sin(2*x)
a

$$a_{1} \cos{\left (x \right )} + a_{2} \cos{\left (2 x \right )} + b_{1} \sin{\left (x \right )} + b_{2} \sin{\left (2 x \right )}$$

We want to square it and get a truncated Fourier series again.

a=(a**2).rewrite(exp).expand().rewrite(cos).expand()
a

$$\frac{a_{1}^{2}}{2} \cos{\left (2 x \right )} + \frac{a_{1}^{2}}{2} + a_{1} a_{2} \cos{\left (x \right )} + a_{1} a_{2} \cos{\left (3 x \right )} + a_{1} b_{1} \sin{\left (2 x \right )} + a_{1} b_{2} \sin{\left (x \right )} + a_{1} b_{2} \sin{\left (3 x \right )} + \frac{a_{2}^{2}}{2} \cos{\left (4 x \right )} + \frac{a_{2}^{2}}{2} - a_{2} b_{1} \sin{\left (x \right )} + a_{2} b_{1} \sin{\left (3 x \right )} + a_{2} b_{2} \sin{\left (4 x \right )} - \frac{b_{1}^{2}}{2} \cos{\left (2 x \right )} + \frac{b_{1}^{2}}{2} + b_{1} b_{2} \cos{\left (x \right )} - b_{1} b_{2} \cos{\left (3 x \right )} - \frac{b_{2}^{2}}{2} \cos{\left (4 x \right )} + \frac{b_{2}^{2}}{2}$$

a.collect([cos(x),cos(2*x),cos(3*x),sin(x),sin(2*x),sin(3*x)])

$$\frac{a_{1}^{2}}{2} + a_{1} b_{1} \sin{\left (2 x \right )} + \frac{a_{2}^{2}}{2} \cos{\left (4 x \right )} + \frac{a_{2}^{2}}{2} + a_{2} b_{2} \sin{\left (4 x \right )} + \frac{b_{1}^{2}}{2} - \frac{b_{2}^{2}}{2} \cos{\left (4 x \right )} + \frac{b_{2}^{2}}{2} + \left(\frac{a_{1}^{2}}{2} - \frac{b_{1}^{2}}{2}\right) \cos{\left (2 x \right )} + \left(a_{1} a_{2} - b_{1} b_{2}\right) \cos{\left (3 x \right )} + \left(a_{1} a_{2} + b_{1} b_{2}\right) \cos{\left (x \right )} + \left(a_{1} b_{2} - a_{2} b_{1}\right) \sin{\left (x \right )} + \left(a_{1} b_{2} + a_{2} b_{1}\right) \sin{\left (3 x \right )}$$

The function expand_log transforms logarithms of products and powers (of positive quantities) into sums of logarithms; logcombine performs the inverse transformation.

a=expand_log(log(p*q**2))
a

$$\log{\left (p \right )} + 2 \log{\left (q \right )}$$

logcombine(a)

$$\log{\left (p q^{2} \right )}$$

The function expand_power_exp rewrites powers whose exponents are sums via products of powers.

expand_power_exp(x**(p+q))

$$x^{p} x^{q}$$

The function expand_power_base rewrites powers whose bases are products via products of powers.

expand_power_base((x*y)**n)

$$x^{n} y^{n}$$

The function powsimp performs the inverse transformations.

powsimp(exp(x)*exp(2*y)),powsimp(x**n*y**n)

$$\left ( e^{x + 2 y}, \quad \left(x y\right)^{n}\right )$$

New symbolic functions can be introduced. They may have an arbitrary numbers of arguments.

f=Function('f')
f(x)+f(x,y)

$$f{\left (x \right )} + f{\left (x,y \right )}$$

Expression structure

Internally expressions are are trees. The function srepr returns a string representing this tree.

srepr(x+1)

"Add(Symbol('x'), Integer(1))"

srepr(x-1)

"Add(Symbol('x'), Integer(-1))"

srepr(x-y)

"Add(Symbol('x'), Mul(Integer(-1), Symbol('y')))"

srepr(2*x*y/3)

"Mul(Rational(2, 3), Symbol('x'), Symbol('y'))"

srepr(x/y)

"Mul(Symbol('x'), Pow(Symbol('y'), Integer(-1)))"

One may use the functions Add, Mul, Pow, etc. instead of the binary operations +, *, **, etc.

Mul(x,Pow(y,-1))

$$\frac{x}{y}$$

srepr(f(x,y))

"Function('f')(Symbol('x'), Symbol('y'))"

The attribute func is the top-level function of an expression, and args is the list of its agruments.

a=2*x*y**2
a.func

sympy.core.mul.Mul

a.args

$$\left ( 2, \quad x, \quad y^{2}\right )$$

for i in a.args:
    print(i)

2
x
y**2

The function subs substitutes an expression for a symbol.

a.subs(y,2)

$$8 x$$

It can perform substitutions for several symbols. To this end, one calls it with a list of tuples or a dictionary.

a.subs([(x,pi),(y,2)])

$$8 \pi$$

a.subs({x:pi,y:2})

$$8 \pi$$

It can substitute not only for a symbol but also for a subexpression - a function with arguments.

a=f(x)+f(y)
a.subs(f(y),1)

$$f{\left (x \right )} + 1$$

(2*x*y*z).subs(x*y,z)

$$2 z^{2}$$

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

$$x^{3} + x + y^{2} + y$$

Substitutions are performed sequentially. In this case, first $x$ is replaced by $y$ producing $y^3+y^2$; then $y$ is replaced by $x$ in this result.

a=x**2+y**3
a.subs([(x,y),(y,x)])

$$x^{3} + x^{2}$$

Interchanging these substitutions leads to a different result.

a.subs([(y,x),(x,y)])

$$y^{3} + y^{2}$$

But if one calls subs with the keyword parameter simultaneous=True, all substitutions are preformed simultaneously. In this way one can, e.g., interchange $x$ and $y$.

a.subs([(x,y),(y,x)],simultaneous=True)

$$x^{3} + y^{2}$$

A function can be replaced by another function.

g=Function('g')
a=f(x)+f(y)
a.subs(f,g)

$$g{\left (x \right )} + g{\left (y \right )}$$

The method replace searches for subexpressions matching a pattern (with wildcards) and replaces them by a given expression.

a=Wild('a')
(f(x)+f(x+y)).replace(f(a),a**2)

$$x^{2} + \left(x + y\right)^{2}$$

(f(x,x)+f(x,y)).replace(f(a,a),a**2)

$$x^{2} + f{\left (x,y \right )}$$

a=x**2+y**2
a.replace(x,x+1)

$$y^{2} + \left(x + 1\right)^{2}$$

Only a complete subtree can match a pattern, not a subset of factors in a product or a smaller power in a larger one.

a=2*x*y*z
a.replace(x*y,z)

$$2 x y z$$

(x+x**2+x**3+x**4).replace(x**2,y)

$$x^{4} + x^{3} + x + y$$

Solving equations

a,b,c,d,e,f=symbols('a b c d e f')

An equation is represented by the function Eq with two arguments. The function solve returns a list of solutions.

solve(Eq(a*x,b),x)

$$\left [ \frac{b}{a}\right ]$$

Instead of equations, one may pass just expressions to solve; they mean equations <expression>=0.

solve(a*x+b,x)

$$\left [ - \frac{b}{a}\right ]$$

A square equation has 2 solutions.

solve(a*x**2+b*x+c,x)

$$\left [ \frac{1}{2 a} \left(- b + \sqrt{- 4 a c + b^{2}}\right), \quad - \frac{1}{2 a} \left(b + \sqrt{- 4 a c + b^{2}}\right)\right ]$$

A system of linear equations.

solve([a*x+b*y-e,c*x+d*y-f],[x,y])

$$\left \{ x : \frac{- b f + d e}{a d - b c}, \quad y : \frac{a f - c e}{a d - b c}\right \}$$

The function roots returns roots of a polynomial together with their multiplicities.

roots(x**3-3*x+2,x)

$$\left \{ -2 : 1, \quad 1 : 2\right \}$$

The function solve_poly_system solves systems of polynomial equations by constructing their Gröbner bases.

p1=x**2+y**2-1
p2=4*x*y-1
solve_poly_system([p1,p2],x,y)

$$\left [ \left ( 4 \left(-1 - \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}}\right) \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}} \left(- \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}} + 1\right), \quad - \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}}\right ), \quad \left ( - 4 \left(-1 + \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}}\right) \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}} \left(\sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}} + 1\right), \quad \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}}\right ), \quad \left ( 4 \left(-1 - \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}}\right) \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}} \left(- \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}} + 1\right), \quad - \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}}\right ), \quad \left ( - 4 \left(-1 + \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}}\right) \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}} \left(\sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}} + 1\right), \quad \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}}\right )\right ]$$

Series

exp(x).series(x,0,5)

$$1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \mathcal{O}\left(x^{5}\right)$$

A series can start from a negative power.

cot(x).series(x,n=5)

$$\frac{1}{x} - \frac{x}{3} - \frac{x^{3}}{45} + \mathcal{O}\left(x^{5}\right)$$

And even run over half-integer powers.

sqrt(x*(1-x)).series(x,n=5)

$$\sqrt{x} - \frac{x^{\frac{3}{2}}}{2} - \frac{x^{\frac{5}{2}}}{8} - \frac{x^{\frac{7}{2}}}{16} - \frac{5 x^{\frac{9}{2}}}{128} + \mathcal{O}\left(x^{5}\right)$$

log(gamma(1+x)).series(x,n=6).rewrite(zeta)

$$- \gamma x + \frac{\pi^{2} x^{2}}{12} - \frac{x^{3} \zeta\left(3\right)}{3} + \frac{\pi^{4} x^{4}}{360} - \frac{x^{5} \zeta\left(5\right)}{5} + \mathcal{O}\left(x^{6}\right)$$

Let's prepare 3 series.

sinx=series(sin(x),x,0,8)
sinx

$$x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \frac{x^{7}}{5040} + \mathcal{O}\left(x^{8}\right)$$

cosx=series(cos(x),x,n=8)
cosx

$$1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} - \frac{x^{6}}{720} + \mathcal{O}\left(x^{8}\right)$$

tanx=series(tan(x),x,n=8)
tanx

$$x + \frac{x^{3}}{3} + \frac{2 x^{5}}{15} + \frac{17 x^{7}}{315} + \mathcal{O}\left(x^{8}\right)$$

Products and ratios of series are not calculated automatically, the function series should be applied to them.

series(tanx*cosx,n=8)

$$x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \frac{x^{7}}{5040} + \mathcal{O}\left(x^{8}\right)$$

series(sinx/cosx,n=8)

$$x + \frac{x^{3}}{3} + \frac{2 x^{5}}{15} + \frac{17 x^{7}}{315} + \mathcal{O}\left(x^{8}\right)$$

And this series should be equal to 1. But since sinx and cosx are known only with a limited accuracy, we obtain 1 with the same accuracy.

series(sinx**2+cosx**2,n=8)

$$1 + \mathcal{O}\left(x^{8}\right)$$

Here the leading terms have canceled, and the result can be obtained only with a lower accuracy.

series((1-cosx)/x**2,n=6)

$$\frac{1}{2} - \frac{x^{2}}{24} + \frac{x^{4}}{720} + \mathcal{O}\left(x^{6}\right)$$

Series can be differentiated and integrated.

diff(sinx,x)

$$1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} - \frac{x^{6}}{720} + \mathcal{O}\left(x^{7}\right)$$

integrate(cosx,x)

$$x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \frac{x^{7}}{5040} + \mathcal{O}\left(x^{9}\right)$$

A series (starting from a small term) can be substituted for an expansion variable in another series. Here are $\sin(\tan(x))$ and $\tan(\sin(x))$.

st=series(sinx.subs(x,tanx),n=8)
st

$$x + \frac{x^{3}}{6} - \frac{x^{5}}{40} - \frac{55 x^{7}}{1008} + \mathcal{O}\left(x^{8}\right)$$

ts=series(tanx.subs(x,sinx),n=8)
ts

$$x + \frac{x^{3}}{6} - \frac{x^{5}}{40} - \frac{107 x^{7}}{5040} + \mathcal{O}\left(x^{8}\right)$$

series(ts-st,n=8)

$$\frac{x^{7}}{30} + \mathcal{O}\left(x^{8}\right)$$

It is not possible to substitute a numerical value for the expansion variable in a series (and hence to plot it). To this end one has to remove the $\mathcal{O}$ term first, transforming a series into a polynomial.

sinx.removeO()

$$- \frac{x^{7}}{5040} + \frac{x^{5}}{120} - \frac{x^{3}}{6} + x$$

Derivatives

a=x*sin(x+y)
diff(a,x)

$$x \cos{\left (x + y \right )} + \sin{\left (x + y \right )}$$

diff(a,y)

$$x \cos{\left (x + y \right )}$$

The second derivative in $x$ and the first one in $y$.

diff(a,x,2,y)

$$- x \cos{\left (x + y \right )} + 2 \sin{\left (x + y \right )}$$

Expressions with undefined functions can be differentiated.

a=x*f(x**2)
b=diff(a,x)
b

$$2 x^{2} \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x^{2} }} + f{\left (x^{2} \right )}$$

What's this?

print(b)

2*x**2*Subs(Derivative(f(_xi_1), _xi_1), (_xi_1,), (x**2,)) + f(x**2)

The function Derivative represents an unevaluated derivative. It can be evaluated by the method doit.

a=Derivative(sin(x),x)
Eq(a,a.doit())

$$\frac{d}{d x} \sin{\left (x \right )} = \cos{\left (x \right )}$$

Integrals

integrate(1/(x*(x**2-2)**2),x)

$$\frac{1}{4} \log{\left (x \right )} - \frac{1}{8} \log{\left (x^{2} - 2 \right )} - \frac{1}{4 x^{2} - 8}$$

integrate(1/(exp(x)+1),x)

$$x - \log{\left (e^{x} + 1 \right )}$$

integrate(log(x),x)

$$x \log{\left (x \right )} - x$$

integrate(x*sin(x),x)

$$- x \cos{\left (x \right )} + \sin{\left (x \right )}$$

integrate(x*exp(-x**2),x)

$$- \frac{e^{- x^{2}}}{2}$$

a=integrate(x**x,x)
a

$$\int x^{x}\, dx$$

This is an unevaluated integral.

print(a)

Integral(x**x, x)

a=Integral(sin(x),x)
Eq(a,a.doit())

$$\int \sin{\left (x \right )}\, dx = - \cos{\left (x \right )}$$

Definite integrals.

integrate(sin(x),(x,0,pi))

$$2$$

oo means $\infty$.

integrate(exp(-x**2),(x,0,oo))

$$\frac{\sqrt{\pi}}{2}$$

integrate(log(x)/(1-x),(x,0,1))

$$- \frac{\pi^{2}}{6}$$

Summing series

summation(1/n**2,(n,1,oo))

$$\frac{\pi^{2}}{6}$$

summation((-1)**n/n**2,(n,1,oo))

$$- \frac{\pi^{2}}{12}$$

summation(1/n**4,(n,1,oo))

$$\frac{\pi^{4}}{90}$$

An unevaluated sum is denoted Sum.

a=Sum(x**n/factorial(n),(n,0,oo))
Eq(a,a.doit())

$$\sum_{n=0}^{\infty} \frac{x^{n}}{n!} = e^{x}$$

Limits

limit((tan(sin(x))-sin(tan(x)))/x**7,x,0)

$$\frac{1}{30}$$

This limit is easy: just expand the numerator and the denominator into series. Things become more difficult if $x=0$ is an essential singularity. Let's calculate one-sided limits.

limit((tan(sin(x))-sin(tan(x)))/(x**7+exp(-1/x)),x,0,'+')

$$\frac{1}{30}$$

limit((tan(sin(x))-sin(tan(x)))/(x**7+exp(-1/x)),x,0,'-')

$$0$$

Differential equations

t=Symbol('t')
x=Function('x')
p=Function('p')

First order.

dsolve(diff(x(t),t)+x(t),x(t))

$$x{\left (t \right )} = C_{1} e^{- t}$$

Second order.

dsolve(diff(x(t),t,2)+x(t),x(t))

$$x{\left (t \right )} = C_{1} \sin{\left (t \right )} + C_{2} \cos{\left (t \right )}$$

A system of first-order equations.

dsolve((diff(x(t),t)-p(t),diff(p(t),t)+x(t)))

$$\left [ x{\left (t \right )} = C_{1} \sin{\left (t \right )} + C_{2} \cos{\left (t \right )}, \quad p{\left (t \right )} = C_{1} \cos{\left (t \right )} - C_{2} \sin{\left (t \right )}\right ]$$

Linear algebra

a,b,c,d,e,f=symbols('a b c d e f')

A matrix can be constructed from a list of lists.

M=Matrix([[a,b,c],[d,e,f]])
M

$$\left[\begin{matrix}a & b & c\\d & e & f\end{matrix}\right]$$

M.shape

$$\left ( 2, \quad 3\right )$$

A row matrix.

Matrix([[1,2,3]])

$$\left[\begin{matrix}1 & 2 & 3\end{matrix}\right]$$

A column matrix.

Matrix([1,2,3])

$$\left[\begin{matrix}1\\2\\3\end{matrix}\right]$$

A matrix can be constructed from a function.

def g(i,j):
    return Rational(1,i+j+1)
Matrix(3,3,g)

$$\left[\begin{matrix}1 & \frac{1}{2} & \frac{1}{3}\\\frac{1}{2} & \frac{1}{3} & \frac{1}{4}\\\frac{1}{3} & \frac{1}{4} & \frac{1}{5}\end{matrix}\right]$$

Or from an undefined function.

g=Function('g')
M=Matrix(3,3,g)
M

$$\left[\begin{matrix}g{\left (0,0 \right )} & g{\left (0,1 \right )} & g{\left (0,2 \right )}\\g{\left (1,0 \right )} & g{\left (1,1 \right )} & g{\left (1,2 \right )}\\g{\left (2,0 \right )} & g{\left (2,1 \right )} & g{\left (2,2 \right )}\end{matrix}\right]$$

M[1,2]

$$g{\left (1,2 \right )}$$

M[1,2]=0
M

$$\left[\begin{matrix}g{\left (0,0 \right )} & g{\left (0,1 \right )} & g{\left (0,2 \right )}\\g{\left (1,0 \right )} & g{\left (1,1 \right )} & 0\\g{\left (2,0 \right )} & g{\left (2,1 \right )} & g{\left (2,2 \right )}\end{matrix}\right]$$

M[2,:]

$$\left[\begin{matrix}g{\left (2,0 \right )} & g{\left (2,1 \right )} & g{\left (2,2 \right )}\end{matrix}\right]$$

M[:,1]

$$\left[\begin{matrix}g{\left (0,1 \right )}\\g{\left (1,1 \right )}\\g{\left (2,1 \right )}\end{matrix}\right]$$

M[0:2,1:3]

$$\left[\begin{matrix}g{\left (0,1 \right )} & g{\left (0,2 \right )}\\g{\left (1,1 \right )} & 0\end{matrix}\right]$$

A unit matrix.

eye(3)

$$\left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$

A zero matrix.

zeros(3)

$$\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]$$

zeros(2,3)

$$\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]$$

A diagonal matrix.

diag(1,2,3)

$$\left[\begin{matrix}1 & 0 & 0\\0 & 2 & 0\\0 & 0 & 3\end{matrix}\right]$$

M=Matrix([[a,1],[0,a]])
diag(1,M,2)

$$\left[\begin{matrix}1 & 0 & 0 & 0\\0 & a & 1 & 0\\0 & 0 & a & 0\\0 & 0 & 0 & 2\end{matrix}\right]$$

Operations with matrices.

A=Matrix([[a,b],[c,d]])
B=Matrix([[1,2],[3,4]])
A+B

$$\left[\begin{matrix}a + 1 & b + 2\\c + 3 & d + 4\end{matrix}\right]$$

A*B,B*A

$$\left ( \left[\begin{matrix}a + 3 b & 2 a + 4 b\\c + 3 d & 2 c + 4 d\end{matrix}\right], \quad \left[\begin{matrix}a + 2 c & b + 2 d\\3 a + 4 c & 3 b + 4 d\end{matrix}\right]\right )$$

A*B-B*A

$$\left[\begin{matrix}3 b - 2 c & 2 a + 3 b - 2 d\\- 3 a - 3 c + 3 d & - 3 b + 2 c\end{matrix}\right]$$

simplify(A**(-1))

$$\left[\begin{matrix}\frac{d}{a d - b c} & - \frac{b}{a d - b c}\\- \frac{c}{a d - b c} & \frac{a}{a d - b c}\end{matrix}\right]$$

det(A)

$$a d - b c$$

Eigenvalues and eigenvectors

x=Symbol('x',real=True)

M=Matrix([[(1-x)**3*(3+x),4*x*(1-x**2),-2*(1-x**2)*(3-x)],
          [4*x*(1-x**2),-(1+x)**3*(3-x),2*(1-x**2)*(3+x)],
          [-2*(1-x**2)*(3-x),2*(1-x**2)*(3+x),16*x]])
M

$$\left[\begin{matrix}\left(- x + 1\right)^{3} \left(x + 3\right) & 4 x \left(- x^{2} + 1\right) & \left(- x + 3\right) \left(2 x^{2} - 2\right)\\4 x \left(- x^{2} + 1\right) & - \left(- x + 3\right) \left(x + 1\right)^{3} & \left(x + 3\right) \left(- 2 x^{2} + 2\right)\\\left(- x + 3\right) \left(2 x^{2} - 2\right) & \left(x + 3\right) \left(- 2 x^{2} + 2\right) & 16 x\end{matrix}\right]$$

det(M)

$$0$$

This means that this matrix has a null space (this matrix transforms vectors from this subspace into 0). Let's find a basis of this subspace.

v=M.nullspace()
len(v)

$$1$$

It is one-dimensional.

v=simplify(v[0])
v

$$\left[\begin{matrix}- \frac{2}{x - 1}\\\frac{2}{x + 1}\\1\end{matrix}\right]$$

Let's check.

simplify(M*v)

$$\left[\begin{matrix}0\\0\\0\end{matrix}\right]$$

The eigenvalues and their multiplicities.

M.eigenvals()

$$\left \{ 0 : 1, \quad - \left(x^{2} + 3\right)^{2} : 1, \quad \left(x^{2} + 3\right)^{2} : 1\right \}$$

If both eigenvalues and corresponding eigenvectors are needed, the method eigenvects is used. It returns a list of tuples. In each tuple the zeroth element is an eigenvalue, the first one is its multiplicity, and the last one is a list of corresponding basis eigenvectors (their number is the multiplicity).

v=M.eigenvects()
len(v)

$$3$$

for i in range(len(v)):
    v[i][2][0]=simplify(v[i][2][0])
v

$$\left [ \left ( 0, \quad 1, \quad \left [ \left[\begin{matrix}- \frac{2}{x - 1}\\\frac{2}{x + 1}\\1\end{matrix}\right]\right ]\right ), \quad \left ( - \left(x^{2} + 3\right)^{2}, \quad 1, \quad \left [ \left[\begin{matrix}\frac{x}{2} + \frac{1}{2}\\\frac{x + 1}{x - 1}\\1\end{matrix}\right]\right ]\right ), \quad \left ( \left(x^{2} + 3\right)^{2}, \quad 1, \quad \left [ \left[\begin{matrix}\frac{x - 1}{x + 1}\\- \frac{x}{2} + \frac{1}{2}\\1\end{matrix}\right]\right ]\right )\right ]$$

Let's check.

for i in range(len(v)):
    z=M*v[i][2][0]-v[i][0]*v[i][2][0]
    pprint(simplify(z))

⎡0⎤
⎢ ⎥
⎢0⎥
⎢ ⎥
⎣0⎦
⎡0⎤
⎢ ⎥
⎢0⎥
⎢ ⎥
⎣0⎦
⎡0⎤
⎢ ⎥
⎢0⎥
⎢ ⎥
⎣0⎦

Jordan normal form

M=Matrix([[Rational(13,9),-Rational(2,9),Rational(1,3),Rational(4,9),Rational(2,3)],
          [-Rational(2,9),Rational(10,9),Rational(2,15),-Rational(2,9),-Rational(11,15)],
          [Rational(1,5),-Rational(2,5),Rational(41,25),-Rational(2,5),Rational(12,25)],
          [Rational(4,9),-Rational(2,9),Rational(14,15),Rational(13,9),-Rational(2,15)],
          [-Rational(4,15),Rational(8,15),Rational(12,25),Rational(8,15),Rational(34,25)]])
M

$$\left[\begin{matrix}\frac{13}{9} & - \frac{2}{9} & \frac{1}{3} & \frac{4}{9} & \frac{2}{3}\\- \frac{2}{9} & \frac{10}{9} & \frac{2}{15} & - \frac{2}{9} & - \frac{11}{15}\\\frac{1}{5} & - \frac{2}{5} & \frac{41}{25} & - \frac{2}{5} & \frac{12}{25}\\\frac{4}{9} & - \frac{2}{9} & \frac{14}{15} & \frac{13}{9} & - \frac{2}{15}\\- \frac{4}{15} & \frac{8}{15} & \frac{12}{25} & \frac{8}{15} & \frac{34}{25}\end{matrix}\right]$$

The method M.jordan_form() returns a couple of matrices, the transformation matrix $P$ and the Jordan form $J$: $M = P J P^{-1}$.

P,J=M.jordan_form()
J

$$\left[\begin{matrix}1 & 0 & 0 & 0 & 0\\0 & 2 & 1 & 0 & 0\\0 & 0 & 2 & 0 & 0\\0 & 0 & 0 & 1 - i & 0\\0 & 0 & 0 & 0 & 1 + i\end{matrix}\right]$$

P=simplify(P)
P

$$\left[\begin{matrix}-2 & \frac{10}{9} & 0 & \frac{5 i}{12} & - \frac{5 i}{12}\\-2 & - \frac{5}{9} & 0 & - \frac{5 i}{6} & \frac{5 i}{6}\\0 & 0 & \frac{4}{3} & - \frac{3}{4} & - \frac{3}{4}\\1 & \frac{10}{9} & 0 & - \frac{5 i}{6} & \frac{5 i}{6}\\0 & 0 & 1 & 1 & 1\end{matrix}\right]$$

Let's check.

Z=P*J*P**(-1)-M
simplify(Z)

$$\left[\begin{matrix}0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\end{matrix}\right]$$

Plots

SymPy uses matplotlib. However, it distributes $x$ points adaptively, not uniformly.

%matplotlib inline

A single function.

plot(sin(x)/x,(x,-10,10))

<sympy.plotting.plot.Plot at 0xaf87e90c>

Several functions.

plot(x,x**2,x**3,(x,0,2))

<sympy.plotting.plot.Plot at 0xb136652c>

Some additional plotting functions can be imported from sympy.plotting.

from sympy.plotting import (plot_parametric,plot_implicit,
                            plot3d,plot3d_parametric_line,
                            plot3d_parametric_surface)

A parametric plot - a Lissajous curve.

t=Symbol('t')
plot_parametric(sin(2*t),cos(3*t),(t,0,2*pi),
                title='Lissajous',xlabel='x',ylabel='y')

<sympy.plotting.plot.Plot at 0xaf709d0c>

An implicit plot - a circle.

plot_implicit(x**2+y**2-1,(x,-1,1),(y,-1,1))

<sympy.plotting.plot.Plot at 0xaf76e1ec>

A surface. If it is not inline but in a separaye window, you can rotate it with your mouse.

plot3d(x*y,(x,-2,2),(y,-2,2))

<sympy.plotting.plot.Plot at 0xaf71ad2c>

Several surfaces.

plot3d(x**2+y**2,x*y,(x,-2,2),(y,-2,2))

<sympy.plotting.plot.Plot at 0xaf5a6f4c>

A parametric space curve - a spiral.

a=0.1
plot3d_parametric_line(cos(t),sin(t),a*t,(t,0,4*pi))

<sympy.plotting.plot.Plot at 0xaf7d962c>

A parametric surface - a torus.

u,v=symbols('u v')
a=0.3
plot3d_parametric_surface((1+a*cos(u))*cos(v),
                          (1+a*cos(u))*sin(v),a*sin(u),
                          (u,0,2*pi),(v,0,2*pi))

<sympy.plotting.plot.Plot at 0xaf664e4c>