import sympy as sp
string = '1/(x-1) + 1/(x+1) + x + 1'

from sympy import init_printing, latex
init_printing()

expr = sp.sympify(string)
sym, = expr.free_symbols
sym

$\displaystyle x$

x = sp.Symbol(sym.name, real=True)

expr = expr.subs(sym,x)
expr

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

frac = sp.cancel(sp.together(expr))
frac

$\displaystyle \frac{x^{3} + x^{2} + x - 1}{x^{2} - 1}$

den = sp.denom(frac)
den

$\displaystyle x^{2} - 1$

poles = sp.solve(den,x)
poles

$\displaystyle \left[ -1, \ 1\right]$

# the full real line is Interval(-oo,oo)
from sympy import oo
sp.Interval(-oo,oo)

$\displaystyle \left(-\infty, \infty\right)$

sp.FiniteSet.fromiter(poles)

$\displaystyle \left\{-1, 1\right\}$

The Domain of Definition

The domain of definition of the function is the set of values of x for which the expression is well-defined, these are values for which the denominator is not zero. So the domain includes everything up to 1 or -1, thats what is meant by the ) instead of ].

domain = sp.Interval(-oo,oo) - sp.FiniteSet.fromiter(poles) 
domain

$\displaystyle \left(-\infty, -1\right) \cup \left(-1, 1\right) \cup \left(1, \infty\right)$

deriv = sp.cancel(sp.diff(frac,x))
deriv

$\displaystyle \frac{x^{4} - 4 x^{2} - 1}{x^{4} - 2 x^{2} + 1}$

extrema = sp.solve(deriv, x)
extrema

$\displaystyle \left[ - \sqrt{2 + \sqrt{5}}, \ \sqrt{2 + \sqrt{5}}\right]$

extrema_values = [frac.subs(x,x0) for x0 in extrema]
extrema_values

$\displaystyle \left[ \frac{- \left(2 + \sqrt{5}\right)^{\frac{3}{2}} - \sqrt{2 + \sqrt{5}} + 1 + \sqrt{5}}{1 + \sqrt{5}}, \ \frac{1 + \sqrt{2 + \sqrt{5}} + \sqrt{5} + \left(2 + \sqrt{5}\right)^{\frac{3}{2}}}{1 + \sqrt{5}}\right]$

m = sp.limit(frac/x,x,oo)
p = sp.limit(frac - m*x,x,oo)
m,p

$\displaystyle \left( 1, \ 1\right)$

m.is_finite

True

def find_asymptote(expr, x):
    """
    Return m,p such that y=m*x+p is an asymptote to the curve y = expr.
    
    If there is no asymptote, return None.
    """
    m = sp.limit(expr/x, x, oo)
    if not m.is_finite:
        return None
    else:
        p = sp.limit(expr-m*x, x, oo)
        return m,p

[find_asymptote(e, x) for e in (frac, x**2, x**2/(x**2+1))]

[(1, 1), None, (0, 1)]

import matplotlib.pyplot as plt
import numpy as np

def numnuts(a,b,c,d,e):
    print(b,c)

numnuts(1,*(2,3),*(4,5))

2 3

def plot_curve(expr, x, x_min, x_max, y_min, y_max):
    """ Plot y = expr(x) over the specified domain"""
    func = sp.lambdify(x, expr, 'numpy')
    xs = np.linspace(x_min, x_max, 200)
    plt.plot(xs, func(xs))
    plt.ylim(y_min, y_max)
    plt.xlim(x_min, x_max)

plot_curve(frac, x, *(-5, 5), *(-10, 10))

[0]+poles+extrema

$\displaystyle \left[ 0, \ -1, \ 1, \ - \sqrt{2 + \sqrt{5}}, \ \sqrt{2 + \sqrt{5}}\right]$

def choose_domain(values):
    values = list(map(float, values))
    low = min(values) - 2
    high = max(values) + 2
    return low,high

values = map(float, [0]+poles+extrema)
min(list(values))

$\displaystyle -2.0581710272714924$

x_min, x_max = choose_domain([0]+poles+extrema)
y_min, y_max = choose_domain([0]+extrema_values+[m*x_min + p, m*x_max + p])
plot_curve(frac, x, *(x_min, x_max), *(y_min, y_max))

plt.vlines(poles, y_min, y_max, 'r', linewidth=2)
plt.plot([x_min, x_max], [m*x_min+p,m*x_max+p], 'r', linewidth=2)
plt.plot(extrema, extrema_values, 'ro')
plt.title("Sketch of $%s \mapsto %s$" % (latex(x), latex(frac)));

Numbers in Expressions

There are several building blocks for expressions the first of which are numbers. There are two kinds of numbers exact expressions and floating point numbers. The latter are represented by the class Float.

• Float;
• Integer;
• Rational.

sp.Integer(3)/sp.Integer(4) # Returns a rational

$\displaystyle \frac{3}{4}$

sp.sympify(3)
sp.S(3)

$\displaystyle 3$

sp.factorint(42)

$\displaystyle \left\{ 2 : 1, \ 3 : 1, \ 7 : 1\right\}$

z = sp.Symbol('z') # a complex variable
x = sp.Symbol('x',real=True) # a real variable
a = sp.Symbol('a', positive=True) # positive and therefore real
c = sp.Symbol('c', constant=True)
type(z)

sympy.core.symbol.Symbol

sp.symbols('x:5')

$\displaystyle \left( x_{0}, \ x_{1}, \ x_{2}, \ x_{3}, \ x_{4}\right)$

Applying a function to an exact number expression always results in another expression. This results in an unevaluated expression.

f = sp.Function('f')
f(0)

$\displaystyle f{\left(0 \right)}$

g = sp.Lambda(x,1-x**2)
g(2)

$\displaystyle -3$

The structure of expressions

Sympy objects are all considered atoms. These are integers, rationals, constants, symbols. All expressions are built on top of these simpler objects. The obey the following invariant expr == expr.class(*expr.args), which means that all their properties derive solely from their type and the contents of their .args attribute. Also since the .args are also sympy objects, the whole expression has a tree structure, which SymPy classes as nodes and atoms as the leaves.

expr = sp.exp(-x)*sp.log(1-x)*2
sp.srepr(expr)

"Mul(Integer(2), exp(Mul(Integer(-1), Symbol('x', real=True))), log(Add(Integer(1), Mul(Integer(-1), Symbol('x', real=True)))))"

Obtain an aggregate view of the contents of an expression with the .atoms() method. By default it returns the set of atomic objects contained in the expressions.

expr.atoms(sp.Integer)

$\displaystyle \left\{-1, 1, 2\right\}$

Substitution

Return a new expression when you do a substitution.

x,y = sp.symbols('x, y', real=True)
(sp.cos(x)*sp.exp(y)).subs({x:sp.pi,y:2})

$\displaystyle - e^{2}$

Checking for Mathematical Equality

To check whether two expressions are mathematically equal, you need to check whether their differences reduce to zero. To do this, apply a canonicalization function to the difference, and test whether its equal to zero. simplify() is generally a useful choice for this but expand() can also be used.

sp.simplify(x*(x+1) - (x**2+x)) == 0

True

Use N() to Evaluate Symbolic Expressions

1. Create a symbolic expression.
2. Do the computation symbollically.
3. Evaluate the result using N() with the subs option.

expr = sp.exp(x)/(x+1)
deriv = sp.diff(expr,x)
sp.N(deriv, subs={x:1.2345})

$\displaystyle 0.849702535502635$

Use lambdify()

To convert symbolic expressions obtained with SymPy into a more efficient form for numerical evaluation.

f = sp.lambdify([x,a],a+x**2,"numpy")
arr = np.random.randn(1000)
result = f(arr, 0.4)
result[:3]

array([0.50639367, 2.0752103 , 0.5414924 ])

Calculus

f = sp.Function('f')
sp.diff(f(x**2),x)

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

sp.limit(sp.exp(-x),x,oo)

$\displaystyle 0$

With a direction to approach the target limit.

sp.limit(1/x,x,0,"-")

$\displaystyle -\infty$

Integrals

Can find most integrals of logarithmic and exponential functions expressible with special functions. Can compute both definite integrals and antiderivatives.

For the anti-derivative of $sin(x)$

sp.integrate(sp.sin(x),x)

$\displaystyle - \cos{\left(x \right)}$

For the definite integral of $sin(x)$ from 0 to pi

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

$\displaystyle 2$

Unevaluated symbolic integrals and antiderivatives are represented by the Integral class. integrate may return these objects if it cannot evaluate the integral. It is also possible to create Integral objects directly, using the same syntax as integrate(). To evaluate them, call their .doit() method.

sp.Integral(sp.sin(x),x)

$\displaystyle \int \sin{\left(x \right)}\, dx$

sp.Integral(sp.sin(x),(x,0,sp.pi))

$\displaystyle \int\limits_{0}^{\pi} \sin{\left(x \right)}\, dx$

i1 = sp.Integral(sp.sin(x),(x,0,sp.pi))
i1.doit()

$\displaystyle 2$

Taylor Series

A Taylor series approximation is an approximation of a function obtained by truncating its Taylor series.

sp.series(sp.cos(x),x)

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

Solving Equations

The main function for solving equations is solve(). This solves by assuming the expression is equal to zero.

sp.solve(x**2-1,x)

$\displaystyle \left[ -1, \ 1\right]$

Solving for a specific variable to break up the expression can be very helpful.

sp.solve(x**2 - sp.exp(a),x)

$\displaystyle \left[ - e^{\frac{a}{2}}, \ e^{\frac{a}{2}}\right]$

To solve a system of equations, pass a list of expressions to solve(): each one will be interpreted, as in the univariate case as an expression of the form expr=0. The result can be returned in one of two forms, depending on the mathematical structure of the input: either as a list of tuples, where each tuple contains the values for the variables in the order given to solve, or a single dictionary for use in subs(), mapping variables to their values.

sp.solve([-sp.exp(x**2)-y,y-3],x,y)

$\displaystyle \left[ \left( - \sqrt{\log{\left(3 \right)} + i \pi}, \ 3\right), \ \left( \sqrt{\log{\left(3 \right)} + i \pi}, \ 3\right)\right]$

Use the dict=True option to always get a dictionary result, that can be used in subs.

sp.solve([-sp.exp(x**2)-y,y-2],x,y,dict=True)

$\displaystyle \left[ \left\{ x : - \sqrt{\log{\left(2 \right)} + i \pi}, \ y : 2\right\}, \ \left\{ x : \sqrt{\log{\left(2 \right)} + i \pi}, \ y : 2\right\}\right]$