

# SymPy

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


```python
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.


```python
x=Symbol('x')
```


```python
a=x**2-1
a
```

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


```python
y,z=symbols('y z')
```

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


```python
a.subs(x,y+1)
```

## Polynomials and rational functions

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


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


```python
a=expand(a)
a
```

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


```python
degree(a,x)
```

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


```python
collect(a,x)
```

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.


```python
a=factor(a)
a
```

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


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


```python
cancel(a)
```

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


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


```python
together(a)
```

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.


```python
simplify(a)
```

Partial fraction decomposition with respect to $x$.


```python
apart(a,x)
```

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


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

And how much is it numerically?


```python
a.n()
```

## Elementary functions

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


```python
sin(-x)
```


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

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


```python
pi.n(200)
```

`E` is the base of natural logarithms.


```python
log(1),log(E)
```


```python
exp(log(x)),log(exp(x))
```

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


```python
sqrt(0)
```


```python
sqrt(x)**4,sqrt(x**4)
```

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.


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


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

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


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

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


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


```python
asin(x).rewrite(log)
```

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


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

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


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

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.


```python
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
```

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


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


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

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


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


```python
logcombine(a)
```

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


```python
expand_power_exp(x**(p+q))
```

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


```python
expand_power_base((x*y)**n)
```

The function `powsimp` performs the inverse transformations.


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

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


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

## Expression structure

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


```python
srepr(x+1)
```




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




```python
srepr(x-1)
```




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




```python
srepr(x-y)
```




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




```python
srepr(2*x*y/3)
```




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




```python
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.


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


```python
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.


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




    sympy.core.mul.Mul




```python
a.args
```


```python
for i in a.args:
    print(i)
```

    2
    x
    y**2


The function `subs` substitutes an expression for a symbol.


```python
a.subs(y,2)
```

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


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


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

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


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


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


```python
(x+x**2+x**3+x**4).subs(x**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.


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

Interchanging these substitutions leads to a different result.


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

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$.


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

A function can be replaced by another function.


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

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


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


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


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

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


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


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

## Solving equations


```python
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.


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

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


```python
solve(a*x+b,x)
```

A square equation has 2 solutions.


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

A system of linear equations.


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

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


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

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


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

## Series


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

A series can start from a negative power.


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

And even run over half-integer powers.


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


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

Let's prepare 3 series.


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


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


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

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


```python
series(tanx*cosx,n=8)
```


```python
series(sinx/cosx,n=8)
```

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.


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

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


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

Series can be differentiated and integrated.


```python
diff(sinx,x)
```


```python
integrate(cosx,x)
```

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))$.


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


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


```python
series(ts-st,n=8)
```

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.


```python
sinx.removeO()
```

## Derivatives


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


```python
diff(a,y)
```

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


```python
diff(a,x,2,y)
```

Expressions with undefined functions can be differentiated.


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

What's this?


```python
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`.


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

## Integrals


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


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


```python
integrate(log(x),x)
```


```python
integrate(x*sin(x),x)
```


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


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

This is an unevaluated integral.


```python
print(a)
```

    Integral(x**x, x)



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

Definite integrals.


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

`oo` means $\infty$.


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


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

## Summing series


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


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


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

An unevaluated sum is denoted `Sum`.


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

## Limits


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

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.


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


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

## Differential equations


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

First order.


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

Second order.


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

A system of first-order equations.


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

## Linear algebra


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

A matrix can be constructed from a list of lists.


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


```python
M.shape
```

A row matrix.


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

A column matrix.


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

A matrix can be constructed from a function.


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

Or from an undefined function.


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


```python
M[1,2]
```


```python
M[1,2]=0
M
```


```python
M[2,:]
```


```python
M[:,1]
```


```python
M[0:2,1:3]
```

A unit matrix.


```python
eye(3)
```

A zero matrix.


```python
zeros(3)
```


```python
zeros(2,3)
```

A diagonal matrix.


```python
diag(1,2,3)
```


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

Operations with matrices.


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


```python
A*B,B*A
```


```python
A*B-B*A
```


```python
simplify(A**(-1))
```


```python
det(A)
```

### Eigenvalues and eigenvectors


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


```python
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
```


```python
det(M)
```

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.


```python
v=M.nullspace()
len(v)
```

It is one-dimensional.


```python
v=simplify(v[0])
v
```

Let's check.


```python
simplify(M*v)
```

The eigenvalues and their multiplicities.


```python
M.eigenvals()
```

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).


```python
v=M.eigenvects()
len(v)
```


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

Let's check.


```python
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


```python
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
```

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


```python
P,J=M.jordan_form()
J
```


```python
P=simplify(P)
P
```

Let's check.


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

## Plots

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


```python
%matplotlib inline
```

A single function.


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

Several functions.


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

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


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

A parametric plot - a Lissajous curve.


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

An implicit plot - a circle.


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

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


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

Several surfaces.


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

A parametric space curve - a spiral.


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

A parametric surface - a torus.


```python
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))
```

183af3da10348f4e367d312987ef698f255c7c88
ipynb
Jupyter Notebook
Gallery/sympy.ipynb
RonSheely/jupyter-notebooks
8d9b37503e59675d392f9b7c3a5a71a1c73a5373
MIT
Gallery/sympy.ipynb
RonSheely/jupyter-notebooks
8d9b37503e59675d392f9b7c3a5a71a1c73a5373
MIT
Gallery/sympy.ipynb
RonSheely/jupyter-notebooks
8d9b37503e59675d392f9b7c3a5a71a1c73a5373
MIT
Qwen/Qwen-72B
1. YES
2. YES