# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # `sympy` guide

import sympy
print(sympy.__version__)

# +
# cf. https://docs.sympy.org/latest/index.html
# Symbolic computation is with symbols, not numerical values
print(sympy.sqrt(8))
print(sympy.sqrt(3))

from sympy import symbols
x, y = symbols('x y')
expr = x + 2*y
print(expr)
print(expr + 1)
print(expr - x) # Simplified when possible; in this case to 2*y
print(x * expr)

from sympy import expand, factor
expanded_expr = expand(x*expr)
print(expanded_expr)
print(factor(expanded_expr))
x, t, z, nu = symbols('x t z nu')
# -

# This will make all further examples pretty print with unicode.
sympy.init_printing(use_unicode=True)

from sympy import (diff, sin, exp)
diff(sin(x) * exp(x), x)

# Compute $\int (e^x \sin{x} + e^x \cos{x}) dx$

from sympy import (cos, integrate)
integrate(exp(x) * sin(x) + exp(x)*cos(x), x)

# Compute $\int^{-\infty}_{\infty} \sin{(x^2)} dx$

from sympy import oo
integrate(sin(x**2), (x, -oo, oo))

# Find $\lim_{x \to 0} \frac{\sin{x}}{x}$

from sympy import limit
limit(sin(x) / x, x, 0)

# ### `solve`, Solve Algebraic equations
#
# Solve $x^2 - 2 = 0$

from sympy import solve
solve(x**2 - 2, x)

# ### `dsolve`, Solve differential equations
# Solve $y'' - y = e^t$

from sympy import (dsolve, Eq, Function)
y = Function('y')
dsolve(Eq(y(t).diff(t, t) - y(t), exp(t)), y(t))

# ### `eigenvals()`, Find eigenvalues

from sympy import Matrix
Matrix([[1, 2], [2, 2]]).eigenvals()

# ### `besselj`, Example of Bessel functions, special functions
#
# Rewrite the Bessel function $J_{\nu}(z)$ in terms of the spherical Bessel function $j_{\nu}(z)$

from sympy import besselj, jn
besselj(nu, z).rewrite(jn)

from sympy import latex, Integral, pi
latex(Integral(cos(x)**2, (x, 0, pi)))

# ## [Sympy Gotchas](https://docs.sympy.org/latest/tutorial/gotchas.html)
#
# Name of a Smybol and name of the variable it's assigned need not have anything to do with one another.
# Best practice usually is to assing Symbols to Python variables of the same name. Symbol names can contain characters that aren't allowed in Python variable names, or assign Symbols long names to single letter Python variables.

crazy = symbols('unrelated')
print(crazy + 1)

# +
x = symbols('x')
expr = x + 1
x = 2
print(expr)

x = 'abc'
expr = x + 'def'
print(expr)
x = 'ABC'
print(expr)

x = symbols('x')
expr = x + 1
print(expr.subs(x, 2))
# -

# Changing `x` to `2` had no effect on `expr`. This is because `x=2` changes the Python variable `x` to `2` but has no effect on SymPy Symbol `x`, which was what we used in creating `expr`. This is how python works.

print(x + 1 == 4)
Eq(x + 1, 4)

# `Eq` is used to create **symbolic equalities**.
#
# It's best to check if a == b by taking a - b and simplify it. It can be [theoretical proven](https://en.wikipedia.org/wiki/Richardson%27s_theorem) that's impossible to determine if 2 symbolic expressions are identically equal in general.

from sympy import simplify
a = (x + 1)**2
b = x**2 + 2*x + 1
print(simplify(a - b))
c = x**2 - 2*x + 1
print(simplify(a - c))

# There's also a method called `equals` that tests if 2 expressions are equal by **evaluting them numerically at random points.**

# +
a = cos(x)**2 - sin(x)**2
b = cos(2*x)
print(a.equals(b))

# EY: 20200909
print(simplify(a - b))
# -

from sympy import Integer
print(Integer(1)/Integer(3))
print(1/2) # PYthon 3 floating point division

from sympy import Rational
print(Rational(1, 2))
print(x + Rational(1, 2))

# Further reading for [Gotchas and Pitfalls](https://docs.sympy.org/latest/gotchas.html#gotchas)

# ## [Basic Operations](https://docs.sympy.org/latest/tutorial/basic_operations.html)

# ### Substitution

# +
x, y, z = symbols("x y z")
expr = cos(x) + 1
print(expr.subs(x, y))
print(expr.subs(x, 0))
expr = x**y
print(expr)
expr = expr.subs(y, x**y)
print(expr)
expr = expr.subs(y, x**x)
print(expr)

from sympy import expand_trig
expr = sin(2*x) + cos(2*x)
print(expand_trig(expr))
# An easy way is just replace sin(2x) with 2sin(x)cos(x)
expr.subs(sin(2*x), 2*sin(x)*cos(x))
print(expr)
# -

# 2 important things to note about `subs`:
# 1. returns new expression, SymPy objects are immutable. That means `subs` doesn't modify it in-place
# 2. 

# +
expr = cos(x)
print(expr.subs(x, 0))
print(expr)
print(x)

expr = x**3 + 4*x*y - z
# To perform multiple substitutions at once, pass a list of (old, new) pairs to subs
print(expr.subs([(x, 2), (y, 4), (z, 0)]))

expr = x**4 - 4*x**3 + 4*x**2 - 2*x + 3
# Often useful to combine this list comprehension to do a large set of
# similar replacements all at once.
# For example, say we wanted to replace all instances of x that have
# an even power with y
replacements = [(x**i, y**i) for i in range(5) if i % 2 == 0]
# -

# #### Converting Strings to SymPy Expressions, `sympify`

from sympy import sympify
str_expr = "x**2 + 3*x - 1/2"
expr = sympify(str_expr)
print(expr)
print(expr.subs(x, 2))

# ## `evalf`, evaluate numerical expression

from sympy import sqrt
expr = sqrt(8) # Default 15 digits of precision
print(expr.evalf())
print(pi.evalf(100)) # Compute first 100 digits of pi

# It's more efficient and numerically stable to pass the substitution to `evalf` using `subs` flag to numerically evaluate an expression.

# +
expr = cos(2*x)
print(expr.evalf(subs={x: 2.4}))

# Set desired precision
one = cos(1)**2 + sin(1)**2
print((one - 1).evalf())
print((one - 1).evalf(chop=True))
# -

# ## `lambdify`, convert Sympy expression to expression for numerics
#
# Convert Sympy expression to an expression that can be numerically evaluated using lambdify.
#
# Warning: lambdify uses eval. Don't use it on unsanitized input.

# +
import numpy
a = numpy.arange(10)
expr = sin(x)
f = lambdify(x, expr, "numpy")
print(f(a))

# Use other libraries than numpy, e.g. use standard library math
f = lambdify(x, expr, "math")
print(f(0.1))

def mysin(x):
    """
    My sine. Note that this is only accurate for small x
    """
    return x
f = lambdify(x, expr, {"sin":mysin})
print(f(0.1))
# -



# Go also to [advanced expression manipulation](https://docs.sympy.org/latest/tutorial/manipulation.html#tutorial-manipulation)





StarcoderdataJupyterScriptsDedupFiltered