## Jupyter Snippet NP ch03-code-listing

Jupyter Snippet NP ch03-code-listing

# Chapter 3: Symbolic computing

Robert Johansson

Source code listings for Numerical Python - Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib (ISBN 978-1-484242-45-2).

import sympy

sympy.init_printing()

from sympy import I, pi, oo

x = sympy.Symbol("x")

y = sympy.Symbol("y", real=True)

y.is_real

True

x.is_real is None

True

sympy.Symbol("z", imaginary=True).is_real

False

x = sympy.Symbol("x")

y = sympy.Symbol("y", positive=True)

sympy.sqrt(x ** 2)


$\displaystyle \sqrt{x^{2}}$

sympy.sqrt(y ** 2)


$\displaystyle y$

n1, n2, n3 = sympy.Symbol("n"), sympy.Symbol("n", integer=True), sympy.Symbol("n", odd=True)

sympy.cos(n1 * pi)


$\displaystyle \cos{\left(\pi n \right)}$

sympy.cos(n2 * pi)


$\displaystyle \left(-1\right)^{n}$

sympy.cos(n3 * pi)


$\displaystyle -1$

a, b, c = sympy.symbols("a, b, c", negative=True)

d, e, f = sympy.symbols("d, e, f", positive=True)


## Numbers

i = sympy.Integer(19)

"i = {} [type {}]".format(i, type(i))

"i = 19 [type <class 'sympy.core.numbers.Integer'>]"

i.is_Integer, i.is_real, i.is_odd

(True, True, True)

f = sympy.Float(2.3)

"f = {} [type {}]".format(f, type(f))

"f = 2.30000000000000 [type <class 'sympy.core.numbers.Float'>]"

f.is_Integer, f.is_real, f.is_odd

(False, True, False)

i, f = sympy.sympify(19), sympy.sympify(2.3)

type(i)

sympy.core.numbers.Integer

type(f)

sympy.core.numbers.Float

n = sympy.Symbol("n", integer=True)

n.is_integer, n.is_Integer, n.is_positive, n.is_Symbol

(True, False, None, True)

i = sympy.Integer(19)

i.is_integer, i.is_Integer, i.is_positive, i.is_Symbol

(True, True, True, False)

i ** 50


$\displaystyle 8663234049605954426644038200675212212900743262211018069459689001$

sympy.factorial(100)


$\displaystyle 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000$

"%.25f" % 0.3  # create a string represention with 25 decimals

'0.2999999999999999888977698'

sympy.Float(0.3, 25)


$\displaystyle 0.2999999999999999888977698$

sympy.Float('0.3', 25)


$\displaystyle 0.3$

### Rationals

sympy.Rational(11, 13)


$\displaystyle \frac{11}{13}$

r1 = sympy.Rational(2, 3)

r2 = sympy.Rational(4, 5)

r1 * r2


$\displaystyle \frac{8}{15}$

r1 / r2


$\displaystyle \frac{5}{6}$

### Functions

x, y, z = sympy.symbols("x, y, z")

f = sympy.Function("f")

type(f)

sympy.core.function.UndefinedFunction

f(x)


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

g = sympy.Function("g")(x, y, z)

g


$\displaystyle g{\left(x,y,z \right)}$

g.free_symbols


$\displaystyle \left{x, y, z\right}$

sympy.sin

sin

sympy.sin(x)


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

sympy.sin(pi * 1.5)


$\displaystyle -1$

n = sympy.Symbol("n", integer=True)

sympy.sin(pi * n)


$\displaystyle 0$

h = sympy.Lambda(x, x**2)

h


$\displaystyle \left( x \mapsto x^{2} \right)$

h(5)


$\displaystyle 25$

h(1+x)


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

### Expressions

x = sympy.Symbol("x")

e = 1 + 2 * x**2 + 3 * x**3

e


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

e.args


$\displaystyle \left( 1, \ 2 x^{2}, \ 3 x^{3}\right)$

e.args[1]


$\displaystyle 2 x^{2}$

e.args[1].args[1]


$\displaystyle x^{2}$

e.args[1].args[1].args[0]


$\displaystyle x$

e.args[1].args[1].args[0].args


$\displaystyle \left( \right)$

## Simplification

expr = 2 * (x**2 - x) - x * (x + 1)

expr


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

sympy.simplify(expr)


$\displaystyle x \left(x - 3\right)$

expr.simplify()


$\displaystyle x \left(x - 3\right)$

expr


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

expr = 2 * sympy.cos(x) * sympy.sin(x)

expr


$\displaystyle 2 \sin{\left(x \right)} \cos{\left(x \right)}$

sympy.trigsimp(expr)


$\displaystyle \sin{\left(2 x \right)}$

expr = sympy.exp(x) * sympy.exp(y)

expr


$\displaystyle e^{x} e^{y}$

sympy.powsimp(expr)


$\displaystyle e^{x + y}$

## Expand

expr = (x + 1) * (x + 2)

sympy.expand(expr)


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

sympy.sin(x + y).expand(trig=True)


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

a, b = sympy.symbols("a, b", positive=True)

sympy.log(a * b).expand(log=True)


$\displaystyle \log{\left(a \right)} + \log{\left(b \right)}$

sympy.exp(I*a + b).expand(complex=True)


$\displaystyle i e^{b} \sin{\left(a \right)} + e^{b} \cos{\left(a \right)}$

sympy.expand((a * b)**x, power_exp=True)


$\displaystyle a^{x} b^{x}$

sympy.exp(I*(a-b)*x).expand(power_exp=True)


$\displaystyle e^{i a x} e^{- i b x}$

## Factor

sympy.factor(x**2 - 1)


$\displaystyle \left(x - 1\right) \left(x + 1\right)$

sympy.factor(x * sympy.cos(y) + sympy.sin(z) * x)


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

sympy.logcombine(sympy.log(a) - sympy.log(b))


$\displaystyle \log{\left(\frac{a}{b} \right)}$

expr = x + y + x * y * z

expr.factor()


$\displaystyle x y z + x + y$

expr.collect(x)


$\displaystyle x \left(y z + 1\right) + y$

expr.collect(y)


$\displaystyle x + y \left(x z + 1\right)$

expr = sympy.cos(x + y) + sympy.sin(x - y)

expr.expand(trig=True).collect([sympy.cos(x), sympy.sin(x)]).collect(sympy.cos(y) - sympy.sin(y))


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

### Together, apart, cancel

sympy.apart(1/(x**2 + 3*x + 2), x)


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

sympy.together(1 / (y * x + y) + 1 / (1+x))


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

sympy.cancel(y / (y * x + y))


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

### Substitutions

(x + y).subs(x, y)


$\displaystyle 2 y$

sympy.sin(x * sympy.exp(x)).subs(x, y)


$\displaystyle \sin{\left(y e^{y} \right)}$

sympy.sin(x * z).subs({z: sympy.exp(y), x: y, sympy.sin: sympy.cos})


$\displaystyle \cos{\left(y e^{y} \right)}$

expr = x * y + z**2 *x

values = {x: 1.25, y: 0.4, z: 3.2}

expr.subs(values)


$\displaystyle 13.3$

## Numerical evaluation

sympy.N(1 + pi)


$\displaystyle 4.14159265358979$

sympy.N(pi, 50)


$\displaystyle 3.1415926535897932384626433832795028841971693993751$

(x + 1/pi).evalf(7)


$\displaystyle x + 0.3183099$

expr = sympy.sin(pi * x * sympy.exp(x))

[expr.subs(x, xx).evalf(3) for xx in range(0, 10)]


$\displaystyle \left[ 0, \ 0.774, \ 0.642, \ 0.722, \ 0.944, \ 0.205, \ 0.974, \ 0.977, \ -0.87, \ -0.695\right]$

expr_func = sympy.lambdify(x, expr)

expr_func(1.0)


$\displaystyle 0.773942685266709$

expr_func = sympy.lambdify(x, expr, 'numpy')

import numpy as np

xvalues = np.arange(0, 10)

expr_func(xvalues)

array([ 0.        ,  0.77394269,  0.64198244,  0.72163867,  0.94361635,
0.20523391,  0.97398794,  0.97734066, -0.87034418, -0.69512687])


## Calculus

f = sympy.Function('f')(x)

sympy.diff(f, x)


$\displaystyle \frac{d}{d x} f{\left(x \right)}$

sympy.diff(f, x, x)


$\displaystyle \frac{d^{2}}{d x^{2}} f{\left(x \right)}$

sympy.diff(f, x, 3)


$\displaystyle \frac{d^{3}}{d x^{3}} f{\left(x \right)}$

g = sympy.Function('g')(x, y)

g.diff(x, y)


$\displaystyle \frac{\partial^{2}}{\partial y\partial x} g{\left(x,y \right)}$

g.diff(x, 3, y, 2)         # equivalent to s.diff(g, x, x, x, y, y)


$\displaystyle \frac{\partial^{5}}{\partial y^{2}\partial x^{3}} g{\left(x,y \right)}$

expr = x**4 + x**3 + x**2 + x + 1

expr.diff(x)


$\displaystyle 4 x^{3} + 3 x^{2} + 2 x + 1$

expr.diff(x, x)


$\displaystyle 2 \left(6 x^{2} + 3 x + 1\right)$

expr = (x + 1)**3 * y ** 2 * (z - 1)

expr.diff(x, y, z)


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

expr = sympy.sin(x * y) * sympy.cos(x / 2)

expr.diff(x)


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

expr = sympy.special.polynomials.hermite(x, 0)

expr.diff(x).doit()


$\displaystyle \frac{2^{x} \sqrt{\pi} \operatorname{polygamma}{\left(0,\frac{1}{2} - \frac{x}{2} \right)}}{2 \Gamma\left(\frac{1}{2} - \frac{x}{2}\right)} + \frac{2^{x} \sqrt{\pi} \log{\left(2 \right)}}{\Gamma\left(\frac{1}{2} - \frac{x}{2}\right)}$

d = sympy.Derivative(sympy.exp(sympy.cos(x)), x)

d


$\displaystyle \frac{d}{d x} e^{\cos{\left(x \right)}}$

d.doit()


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

## Integrals

a, b = sympy.symbols("a, b")
x, y = sympy.symbols('x, y')
f = sympy.Function('f')(x)

sympy.integrate(f)


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

sympy.integrate(f, (x, a, b))


$\displaystyle \int\limits_{a}^{b} f{\left(x \right)}, dx$

sympy.integrate(sympy.sin(x))


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

sympy.integrate(sympy.sin(x), (x, a, b))


$\displaystyle \cos{\left(a \right)} - \cos{\left(b \right)}$

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


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

a, b, c = sympy.symbols("a, b, c", positive=True)

sympy.integrate(a * sympy.exp(-((x-b)/c)**2), (x, -oo, oo))


$\displaystyle \sqrt{\pi} a c$

sympy.integrate(sympy.sin(x * sympy.cos(x)))


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

expr = sympy.sin(x*sympy.exp(y))

sympy.integrate(expr, x)


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

expr = (x + y)**2

sympy.integrate(expr, x)


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

sympy.integrate(expr, x, y)


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

sympy.integrate(expr, (x, 0, 1), (y, 0, 1))


$\displaystyle \frac{7}{6}$

## Series

x = sympy.Symbol("x")

f = sympy.Function("f")(x)

sympy.series(f, x)


$\displaystyle f{\left(0 \right)} + x \left. \frac{d}{d x} f{\left(x \right)} \right|{\substack{ x=0 }} + \frac{x^{2} \left. \frac{d^{2}}{d x^{2}} f{\left(x \right)} \right|{\substack{ x=0 }}}{2} + \frac{x^{3} \left. \frac{d^{3}}{d x^{3}} f{\left(x \right)} \right|{\substack{ x=0 }}}{6} + \frac{x^{4} \left. \frac{d^{4}}{d x^{4}} f{\left(x \right)} \right|{\substack{ x=0 }}}{24} + \frac{x^{5} \left. \frac{d^{5}}{d x^{5}} f{\left(x \right)} \right|_{\substack{ x=0 }}}{120} + O\left(x^{6}\right)$

x0 = sympy.Symbol("{x_0}")

f.series(x, x0, n=2)


$\displaystyle f{\left({x_0} \right)} + \left(x - {x_0}\right) \left. \frac{d}{d \xi_{1}} f{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}={x_0} }} + O\left(\left(x - {x_0}\right)^{2}; x\rightarrow {x_0}\right)$

f.series(x, x0, n=2).removeO()


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

sympy.cos(x).series()


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

sympy.sin(x).series()


$\displaystyle x - \frac{x^{3}}{6} + \frac{x^{5}}{120} + O\left(x^{6}\right)$

sympy.exp(x).series()


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

(1/(1+x)).series()


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

expr = sympy.cos(x) / (1 + sympy.sin(x * y))

expr.series(x, n=4)


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

expr.series(y, n=4)


$\displaystyle \cos{\left(x \right)} - x y \cos{\left(x \right)} + x^{2} y^{2} \cos{\left(x \right)} - \frac{5 x^{3} y^{3} \cos{\left(x \right)}}{6} + O\left(y^{4}\right)$

expr.series(y).removeO().series(x).removeO().expand()


$\displaystyle - \frac{61 x^{5} y^{5}}{120} + \frac{5 x^{5} y^{3}}{12} - \frac{x^{5} y}{24} + \frac{2 x^{4} y^{4}}{3} - \frac{x^{4} y^{2}}{2} + \frac{x^{4}}{24} - \frac{5 x^{3} y^{3}}{6} + \frac{x^{3} y}{2} + x^{2} y^{2} - \frac{x^{2}}{2} - x y + 1$

## Limits

sympy.limit(sympy.sin(x) / x, x, 0)


$\displaystyle 1$

f = sympy.Function('f')
x, h = sympy.symbols("x, h")

diff_limit = (f(x + h) - f(x))/h

sympy.limit(diff_limit.subs(f, sympy.cos), h, 0)


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

sympy.limit(diff_limit.subs(f, sympy.sin), h, 0)


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

expr = (x**2 - 3*x) / (2*x - 2)

p = sympy.limit(expr/x, x, oo)

q = sympy.limit(expr - p*x, x, oo)

p, q


$\displaystyle \left( \frac{1}{2}, \ -1\right)$

## Sums and products

n = sympy.symbols("n", integer=True)

x = sympy.Sum(1/(n**2), (n, 1, oo))

x


$\displaystyle \sum_{n=1}^{\infty} \frac{1}{n^{2}}$

x.doit()


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

x = sympy.Product(n, (n, 1, 7))

x


$\displaystyle \prod_{n=1}^{7} n$

x.doit()


$\displaystyle 5040$

x = sympy.Symbol("x")

sympy.Sum((x)**n/(sympy.factorial(n)), (n, 1, oo)).doit().simplify()


$\displaystyle e^{x} - 1$

## Equations

x = sympy.symbols("x")

sympy.solve(x**2 + 2*x - 3)


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

a, b, c = sympy.symbols("a, b, c")

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


$\displaystyle \left[ \frac{- b + \sqrt{- 4 a c + b^{2}}}{2 a}, \ - \frac{b + \sqrt{- 4 a c + b^{2}}}{2 a}\right]$

sympy.solve(sympy.sin(x) - sympy.cos(x), x)


$\displaystyle \left[ - \frac{3 \pi}{4}, \ \frac{\pi}{4}\right]$

sympy.solve(sympy.exp(x) + 2 * x, x)


$\displaystyle \left[ - \operatorname{LambertW}{\left(\frac{1}{2} \right)}\right]$

sympy.solve(x**5 - x**2 + 1, x)


$\displaystyle \left[ \operatorname{CRootOf} {\left(x^{5} - x^{2} + 1, 0\right)}, \ \operatorname{CRootOf} {\left(x^{5} - x^{2} + 1, 1\right)}, \ \operatorname{CRootOf} {\left(x^{5} - x^{2} + 1, 2\right)}, \ \operatorname{CRootOf} {\left(x^{5} - x^{2} + 1, 3\right)}, \ \operatorname{CRootOf} {\left(x^{5} - x^{2} + 1, 4\right)}\right]$

1 #s.solve(s.tan(x) - x, x)


$\displaystyle 1$

eq1 = x + 2 * y - 1
eq2 = x - y + 1

sympy.solve([eq1, eq2], [x, y], dict=True)


$\displaystyle \left[ \left{ x : - \frac{1}{3}, \ y : \frac{2}{3}\right}\right]$

eq1 = x**2 - y
eq2 = y**2 - x

sols = sympy.solve([eq1, eq2], [x, y], dict=True)

sols


$\displaystyle \left[ \left{ x : 0, \ y : 0\right}, \ \left{ x : 1, \ y : 1\right}, \ \left{ x : \left(- \frac{1}{2} - \frac{\sqrt{3} i}{2}\right)^{2}, \ y : - \frac{1}{2} - \frac{\sqrt{3} i}{2}\right}, \ \left{ x : \left(- \frac{1}{2} + \frac{\sqrt{3} i}{2}\right)^{2}, \ y : - \frac{1}{2} + \frac{\sqrt{3} i}{2}\right}\right]$

[eq1.subs(sol).simplify() == 0 and eq2.subs(sol).simplify() == 0 for sol in sols]

[True, True, True, True]


## Linear algebra

sympy.Matrix([1,2])


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

sympy.Matrix([[1,2]])


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

sympy.Matrix([[1, 2], [3, 4]])


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

sympy.Matrix(3, 4, lambda m,n: 10 * m + n)


$\displaystyle \left[\begin{matrix}0 & 1 & 2 & 3\10 & 11 & 12 & 13\20 & 21 & 22 & 23\end{matrix}\right]$

a, b, c, d = sympy.symbols("a, b, c, d")

M = sympy.Matrix([[a, b], [c, d]])

M


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

M * M


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

x = sympy.Matrix(sympy.symbols("x_1, x_2"))

M * x


$\displaystyle \left[\begin{matrix}a x_{1} + b x_{2}\c x_{1} + d x_{2}\end{matrix}\right]$

p, q = sympy.symbols("p, q")

M = sympy.Matrix([[1, p], [q, 1]])

M


$\displaystyle \left[\begin{matrix}1 & p\q & 1\end{matrix}\right]$

b = sympy.Matrix(sympy.symbols("b_1, b_2"))

b


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

x = M.solve(b)
x


$\displaystyle \left[\begin{matrix}\frac{b_{1} \left(- p q + 1\right) - p \left(- b_{1} q + b_{2}\right)}{- p q + 1}\\frac{- b_{1} q + b_{2}}{- p q + 1}\end{matrix}\right]$

x = M.LUsolve(b)

x


$\displaystyle \left[\begin{matrix}b_{1} - \frac{p \left(- b_{1} q + b_{2}\right)}{- p q + 1}\\frac{- b_{1} q + b_{2}}{- p q + 1}\end{matrix}\right]$

x = M.inv() * b

x


$\displaystyle \left[\begin{matrix}\frac{b_{1}}{- p q + 1} - \frac{b_{2} p}{- p q + 1}\- \frac{b_{1} q}{- p q + 1} + \frac{b_{2}}{- p q + 1}\end{matrix}\right]$

## Versions

%reload_ext version_information
%version_information sympy, numpy