#!/usr/bin/env python
# coding: utf-8

# ![python](python_integral.png)
# 
# # SymPy
# 
# `SymPy` is a package for symbolic calculations in python, similar to *Mathematica*. It works with expressions containing symbols.

# In[1]:


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.

# In[2]:


x=Symbol('x')


# In[3]:


a=x**2-1
a


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

# In[4]:


y,z=symbols('y z')


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

# In[5]:


a.subs(x,y+1)


# ## Polynomials and rational functions
# 
# `SymPy` does not expand brackets automatically. The function `expand` is used for this.

# In[6]:


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


# In[7]:


a=expand(a)
a


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

# In[8]:


degree(a,x)


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

# In[9]:


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.

# In[10]:


a=factor(a)
a


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

# In[11]:


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


# In[12]:


cancel(a)


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

# In[13]:


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


# In[14]:


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.

# In[15]:


simplify(a)


# Partial fraction decomposition with respect to $x$.

# In[16]:


apart(a,x)


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

# In[17]:


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


# And how much is it numerically?

# In[18]:


a.n()


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

# In[19]:


sin(-x)


# In[20]:


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.

# In[21]:


pi.n(100)


# `E` is the base of natural logarithms.

# In[22]:


log(1),log(E)


# In[23]:


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


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

# In[24]:


sqrt(0)


# In[25]:


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.

# In[26]:


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


# In[27]:


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


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

# In[28]:


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


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

# In[29]:


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


# In[30]:


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.

# In[31]:


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


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

# In[32]:


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.

# In[33]:


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.

# In[34]:


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


# In[35]:


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.

# In[36]:


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


# In[37]:


logcombine(a)


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

# In[38]:


expand_power_exp(x**(p+q))


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

# In[39]:


expand_power_base((x*y)**n)


# The function `powsimp` performs the inverse transformations.

# In[40]:


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.

# In[41]:


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


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

# In[42]:


srepr(x+1)


# In[43]:


srepr(x-1)


# In[44]:


srepr(x-y)


# In[45]:


srepr(2*x*y/3)


# In[46]:


srepr(x/y)


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

# In[47]:


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


# In[48]:


srepr(f(x,y))


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

# In[49]:


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


# In[50]:


a.args


# In[51]:


for i in a.args:
    print(i)


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

# In[52]:


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.

# In[53]:


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


# In[54]:


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


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

# In[55]:


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


# In[56]:


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


# In[57]:


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

# In[58]:


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


# Interchanging these substitutions leads to a different result.

# In[59]:


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

# In[60]:


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


# A function can be replaced by another function.

# In[61]:


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.

# In[62]:


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


# In[63]:


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


# In[64]:


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.

# In[65]:


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


# In[66]:


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


# ## Solving equations

# In[67]:


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.

# In[68]:


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


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

# In[69]:


solve(a*x+b,x)


# A square equation has 2 solutions.

# In[70]:


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


# A system of linear equations.

# In[71]:


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.

# In[72]:


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


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

# In[73]:


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


# ## Series

# In[74]:


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


# A series can start from a negative power.

# In[75]:


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


# And even run over half-integer powers.

# In[76]:


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


# In[77]:


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


# Let's prepare 3 series.

# In[78]:


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


# In[79]:


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


# In[80]:


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.

# In[81]:


series(tanx*cosx,n=8)


# In[82]:


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.

# In[83]:


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


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

# In[84]:


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


# Series can be differentiated and integrated.

# In[85]:


diff(sinx,x)


# In[86]:


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

# In[87]:


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


# In[88]:


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


# In[89]:


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.

# In[90]:


sinx.removeO()


# ## Derivatives

# In[91]:


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


# In[92]:


diff(a,y)


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

# In[93]:


diff(a,x,2,y)


# Expressions with undefined functions can be differentiated.

# In[94]:


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


# What's this?

# In[95]:


print(b)


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

# In[96]:


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


# ## Integrals

# In[97]:


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


# In[98]:


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


# In[99]:


integrate(log(x),x)


# In[100]:


integrate(x*sin(x),x)


# In[101]:


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


# In[102]:


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


# This is an unevaluated integral.

# In[103]:


print(a)


# In[104]:


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


# Definite integrals.

# In[105]:


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


# `oo` means $\infty$.

# In[106]:


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


# In[107]:


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


# ## Summing series

# In[108]:


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


# In[109]:


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


# In[110]:


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


# An unevaluated sum is denoted `Sum`.

# In[111]:


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


# ## Limits

# In[112]:


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.

# In[113]:


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


# In[114]:


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


# ## Differential equations

# In[115]:


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


# First order.

# In[116]:


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


# Second order.

# In[117]:


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


# A system of first-order equations.

# In[118]:


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


# ## Linear algebra

# In[119]:


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


# A matrix can be constructed from a list of lists.

# In[120]:


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


# In[121]:


M.shape


# A row matrix.

# In[122]:


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


# A column matrix.

# In[123]:


Matrix([1,2,3])


# A matrix can be constructed from a function.

# In[124]:


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


# Or from an undefined function.

# In[125]:


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


# In[126]:


M[1,2]


# In[127]:


M[1,2]=0
M


# In[128]:


M[2,:]


# In[129]:


M[:,1]


# In[130]:


M[0:2,1:3]


# A unit matrix.

# In[131]:


eye(3)


# A zero matrix.

# In[132]:


zeros(3)


# In[133]:


zeros(2,3)


# A diagonal matrix.

# In[134]:


diag(1,2,3)


# In[135]:


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


# Operations with matrices.

# In[136]:


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


# In[137]:


A*B,B*A


# In[138]:


A*B-B*A


# In[139]:


simplify(A**(-1))


# In[140]:


det(A)


# ### Eigenvalues and eigenvectors

# In[141]:


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


# In[142]:


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


# In[143]:


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.

# In[144]:


v=M.nullspace()
len(v)


# It is one-dimensional.

# In[145]:


v=simplify(v[0])
v


# Let's check.

# In[146]:


simplify(M*v)


# The eigenvalues and their multiplicities.

# In[147]:


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

# In[148]:


v=M.eigenvects()
len(v)


# In[149]:


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


# Let's check.

# In[150]:


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


# ### Jordan normal form

# In[151]:


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

# In[152]:


P,J=M.jordan_form()
J


# In[153]:


P=simplify(P)
P


# Let's check.

# In[154]:


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


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

# In[155]:


get_ipython().run_line_magic('matplotlib', 'inline')


# A single function.

# In[156]:


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


# Several functions.

# In[157]:


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


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

# In[158]:


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


# A parametric plot - a Lissajous curve.

# In[159]:


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.

# In[160]:


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.

# In[161]:


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


# Several surfaces.

# In[162]:


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


# A parametric space curve - a spiral.

# In[163]:


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


# A parametric surface - a torus.

# In[164]:


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