# Computer Algebra System (CAS)

SymPy is a Python library for symbolic mathematics.
- [tutorial](https://docs.sympy.org/latest/tutorial/index.html)
- [documentation](https://docs.sympy.org/latest/tutorial/index.html)

Symbolic calculations begin by defining symbolic variables:


```python
import numpy as np
from sympy import *
x,y=symbols('x y')
```

Evaluate expressions containing symbolic variables:


```python
x+y
```




$\displaystyle x + y$




```python
(x+y)**2
```




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




```python
expand((x + y)**2)
```




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



Substitute variables in expressions: Substitution replaces all instances of something in an expression with something else. It is done using the `subs` method.


```python
z = x**2 + y**2
```


```python
z
```




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



replace y with a + b in the expression for z


```python
a,b=symbols('a b')
z.subs(y,a+b)
```




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



replace x with a + b in the expression for z


```python
z.subs(x,a+b)
```




$\displaystyle y^{2} + \left(a + b\right)^{2}$



Define and evaluate symbolic functions


```python
f = Function('f')
f(x)+1
```




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




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




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



To treat numbers symbolically rather than numerically, sym must be inserted:


```python
sin(np.pi)
```




$\displaystyle 1.22464679914735 \cdot 10^{-16}$




```python
sin(pi)
```




$\displaystyle 0$



Sometimes simplify needs to be called to simplify an expression:


```python
g = (1 + sqrt(Symbol('5')))/2  #golden ratio
g
```




$\displaystyle \frac{\sqrt{5}}{2} + \frac{1}{2}$




```python
g**2-g-1 #should vanish
```




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




```python
simplify(g**2-g-1)
```




$\displaystyle \frac{5}{4} - \frac{5}{4}$



Conversion between symbolic and numeric quantities:


```python
nsimplify(0.25) #symbolic
```




$\displaystyle \frac{1}{4}$




```python
1/3  #numeric
```




    0.3333333333333333




```python
S(1)/3         #symbolic
```




$\displaystyle \frac{1}{3}$




```python
Rational(1,3)  #symbolic
```




$\displaystyle \frac{1}{3}$




```python
Rational(1,3).evalf()
```




$\displaystyle 0.333333333333333$



**Matrices**: A matrix can contain any symbolic expression.


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




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




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




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




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




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




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




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




```python
a,b,c=symbols('a b c')
Matrix([a,b,c])
```




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




```python
Matrix([[1,2],[3,4]])@Matrix([a,b]) # Python 3.5+
```




$\displaystyle \left[\begin{matrix}a + 2 b\\3 a + 4 b\end{matrix}\right]$




```python
Matrix([[1,2],[3,4]])*Matrix([a,b])
```




$\displaystyle \left[\begin{matrix}a + 2 b\\3 a + 4 b\end{matrix}\right]$



calculating the jacobian


```python
x1,x2,x3 = symbols('x1 x2 x3')
f1,f2,f3 = symbols('f1 f2 f3', cls=Function)
M = Matrix([f1(x1)+f1(x2)+f1(x3),f2(x1)+f2(x2)+f2(x3),f3(x1)+f3(x2)+f3(x3)])
M.jacobian([x1,x2,x3])
```




$\displaystyle \left[\begin{matrix}\frac{d}{d x_{1}} \operatorname{f_{1}}{\left(x_{1} \right)} & \frac{d}{d x_{2}} \operatorname{f_{1}}{\left(x_{2} \right)} & \frac{d}{d x_{3}} \operatorname{f_{1}}{\left(x_{3} \right)}\\\frac{d}{d x_{1}} \operatorname{f_{2}}{\left(x_{1} \right)} & \frac{d}{d x_{2}} \operatorname{f_{2}}{\left(x_{2} \right)} & \frac{d}{d x_{3}} \operatorname{f_{2}}{\left(x_{3} \right)}\\\frac{d}{d x_{1}} \operatorname{f_{3}}{\left(x_{1} \right)} & \frac{d}{d x_{2}} \operatorname{f_{3}}{\left(x_{2} \right)} & \frac{d}{d x_{3}} \operatorname{f_{3}}{\left(x_{3} \right)}\end{matrix}\right]$



**Matrix Symbols**: operation on matrices of symbolic dimension $n\times m$


```python
n, m = symbols('n m', integer=True)
M = MatrixSymbol('M', n, m)
b = MatrixSymbol('b', m, 1)
M@b
```




$\displaystyle M b$




```python
(M@b).shape
```




    (n, 1)




```python
Transpose(M@b)
```




$\displaystyle \left(M b\right)^{T}$




```python
Transpose(M@b).doit()
#
#
# ## Algebraic equations
#
# solves the equation for y
```




$\displaystyle b^{T} M^{T}$




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


```python

```


```python
simplify(1/3)
```




$\displaystyle 0.333333333333333$




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




    [-sqrt(y), sqrt(y)]



### Linear systems

Example: $x +2y = u$, $4x + 5y = v$


```python
x, y, u, v = symbols('x y u v')
system = Matrix(( (1, 2, u), (4, 5, v)))
solve_linear_system(system, x, y)
```




    {x: -5*u/3 + 2*v/3, y: 4*u/3 - v/3}




```python
x, y, u, v = symbols('x y u v')
solve((x+2*y-u,4*x+5*y-v), x, y)
```




    {x: -5*u/3 + 2*v/3, y: 4*u/3 - v/3}




```python
x_1, x_2, u, v = symbols('x_1 x_2 u v')
A = Matrix(((1, 2), (4, 5)))
b = Matrix(((u), (v)))
x = Matrix(((x_1), (x_2)))
solve(A*x-b, x_1, x_2)
```




    {x_1: -5*u/3 + 2*v/3, x_2: 4*u/3 - v/3}



### Nonlinear systems

Example: $x^2 + 3y = 0$, $y^2+2x=0$


```python
x, y = symbols('x y')
solve([x**2 + 3*y, y**2 + 2*x])
```




    [{x: 0, y: 0},
     {x: -(12**(1/3)/2 + 2**(2/3)*3**(5/6)*I/2)**2/2,
      y: 12**(1/3)/2 + 2**(2/3)*3**(5/6)*I/2},
     {x: -(-3*2**(2/3)*3**(1/3)/4 - 12**(1/3)/4)**2/2,
      y: -3*2**(2/3)*3**(1/3)/4 - 12**(1/3)/4},
     {x: -(-12**(1/3)/4 + 3*2**(2/3)*3**(1/3)/4 - 2**(2/3)*3**(5/6)*I/2)**2/2,
      y: -12**(1/3)/4 + 3*2**(2/3)*3**(1/3)/4 - 2**(2/3)*3**(5/6)*I/2}]



### Exercise: 

Find all matrices of the form

$$
\left( \begin{array}{ccc} a&b&c\\c&a&b\\b&c&a \end{array} \right)
$$

with $a,b,c\in\mathbf{R}$ whose square is the unit matrix.


```python
a, b, c = symbols('a b c')
M = Matrix(((a, b, c), (c, a, b), (b, c, a)))
solve(M**2-eye(3,3))
```




    [{a: -1, b: 0, c: 0},
     {a: -1/3, b: 2/3, c: 2/3},
     {a: -1/3, b: -1/3 - sqrt(3)*I/3, c: -1/3 + sqrt(3)*I/3},
     {a: -1/3, b: -1/3 + sqrt(3)*I/3, c: -1/3 - sqrt(3)*I/3},
     {a: 1/3, b: -2/3, c: -2/3},
     {a: 1/3, b: 1/3 - sqrt(3)*I/3, c: 1/3 + sqrt(3)*I/3},
     {a: 1/3, b: 1/3 + sqrt(3)*I/3, c: 1/3 - sqrt(3)*I/3},
     {a: 1, b: 0, c: 0}]



### Assumptions on symbolic variables


```python
x = symbols('x')
solve(x**3 + 1)
```




    [-1, 1/2 - sqrt(3)*I/2, 1/2 + sqrt(3)*I/2]




```python
x = symbols('x', real=True)
solve(x**3 + 1)
```




    [-1]




```python
x = symbols('x', positive=True)
solve(x**3 + 1)
```




    []




```python
x.is_real
```




    True




```python
x.is_positive
```




    True



reset assumptions


```python
x = symbols('x')
solve(x**3 + 1)
```




    [-1, 1/2 - sqrt(3)*I/2, 1/2 + sqrt(3)*I/2]



## Symbolic differentiation and integration

Differentiate a symbolic function:


```python
f = Function('f')
Derivative(f(x),x,x)
```




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




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




$\displaystyle 2 x y$




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




$\displaystyle x^{2}$



Second order derivative


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




$\displaystyle 2 y$




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




$\displaystyle 2 y$



### Exercise: 

Determine

$$
\lim_{x\to 0} \frac{2\tan(x) - \tan(2x)}{x(1 - \cos(3x))}
$$

Here we use derivatives to determine limits. This is in particular useful, if we obtain limits that end up in inteterminate or undefined forms such as $\infty / \infty$, $-\infty / \infty$, or $0/0$. The latter is the case here.

#### L'Hopital's Rule:

If $\lim_{x\to 0} f(x)=0$ and $\lim_{x\to 0} g(x)=0$ and $\lim_{x\to 0} \frac{f'(x)}{g'(x)}=L$ then we can say that the limit $\lim_{x\to 0} \frac{f(x)}{g(x)}=L$.


```python
x = symbols('x')
f=(2*tan(x) - tan(2*x))/(x*(1 - cos(3*x)))
f
```




$\displaystyle \frac{2 \tan{\left(x \right)} - \tan{\left(2 x \right)}}{x \left(1 - \cos{\left(3 x \right)}\right)}$




```python
z = 2*tan(x) - tan(2*x)  #numerator
z.subs(x, 0)
```




$\displaystyle 0$




```python
z1 = diff(z)
z1.subs(x,0)
```




$\displaystyle 0$




```python
z2 = diff(z1)
z2.subs(x,0)
```




$\displaystyle 0$




```python
z3 = diff(z2)
z3.subs(x,0)
```




$\displaystyle -12$




```python
n = x*(1 - cos(3*x))  #denominator
n.subs(x, 0)
```




$\displaystyle 0$




```python
n1 = diff(n)
n1.subs(x, 0)
```




$\displaystyle 0$




```python
n2 = diff(n1)
n2.subs(x, 0)
```




$\displaystyle 0$




```python
n3 = diff(n2)
n3.subs(x, 0)
```




$\displaystyle 27$



Applying [L'Hopital's rule](https://www.khanacademy.org/math/ap-calculus-ab/ab-diff-contextual-applications-new/ab-4-7/v/introduction-to-l-hopital-s-rule) three times yields the desired limit:


```python
(z1/n1).subs(x, 0)
```




$\displaystyle \text{NaN}$




```python
(z2/n2).subs(x, 0)
```




$\displaystyle \text{NaN}$




```python
(z3/n3).subs(x, 0)
```




$\displaystyle - \frac{4}{9}$



## Taylor expansion

A [Taylor series](https://mathworld.wolfram.com/TaylorSeries.html) is a series expansion of a function about a point. A one-dimensional Taylor series is an expansion of a real function $f(x)$ about a point $x=a$ is given by

$$f(x)=f(a)+f'(a)(x-a)+\frac{f''(a)}{2!}(x-a)^2+\frac{f'''(a)}{3!}(x-a)^3+\dots+\frac{f^{(n)}(a)}{n!}(x-a)^n+\dots$$

The Taylor series of an unknown function $f(x)$ about a point $a$ up to order $n$ may be found using `series`.


```python
x, h = symbols("x,h")
f = Function("f")
f(x).series(x, x0=h, n=3)
```




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



The Taylor series of the function $\sin(x)$ about a point $x_0=0$ up to order $n=6$:


```python
x = symbols("x")
sin(x).series(x, x0=0, n=6)
```




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



The Taylor series of the function $\sin(x)$ about a point $x_0=0$ up to order $n=10$:


```python
x = symbols("x")
sin(x).series(x, x0=0, n=10)
```




$\displaystyle x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \frac{x^{7}}{5040} + \frac{x^{9}}{362880} + O\left(x^{10}\right)$



The Taylor series of the function $\log(x)$ about a point $x_0=1$ up to order $n=6$:


```python
x = symbols("x")
log(x).series(x, x0=1, n=6)
```




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



If you want a finite (forward) [difference approximation](https://docs.sympy.org/latest/special_topics/finite_diff_derivatives.html) to derivatives, you can for example write


```python
FW = f(x+h).series(x+h, x0=h, n=3)   #TODO
FW.subs(x+h,0)
```




$\displaystyle f{\left(2 h + x \right)}$



## Integration

SymPy implements methods to calculate definite and indefinite [integrals](https://docs.sympy.org/latest/modules/integrals/integrals.html) of expressions.
- `integrate(f,x)` returns the indefinite integral
- `integrate(f,(x,a))` returns the indefinite integral, with result given with $a$ replacing $x$
- `integrate(f,(x,a,b))` returns the definite integral


```python
x = symbols("x")
f = Function("f")
f(x).integrate(x)
```




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




```python
x, a = symbols("x,a")
f = Function("f")
f(x).integrate((x,a))
```




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




```python
x, a, b = symbols("x,a,b")
f = Function("f")
f(x).integrate((x,a,b))
```




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



We start with a quadratic polynomical and take a derivative:


```python
x, a, b, c = symbols("x,a,b,c")
p = a*x**2 + b*x + c
dpdx = diff(p,x)
dpdx
```




$\displaystyle 2 a x + b$



Now we integrate to try and retrieve the function we started with:


```python
p1 = integrate(dpdx,x)
p1
```




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



**Note**: the integration constant is zero by default and thus needs to be added in by hand


```python
p1 = p1 + c
p1
```




$\displaystyle a x^{2} + b x + c$



#### Indefinite, definite and improper integrals

$$f(x)=\int \frac{1}{x^2}$$


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




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




```python
plot(1/x**2, (x, 1, 2))
```


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




$\displaystyle \frac{1}{2}$




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




$\displaystyle 1$



$$\int_0^\inf \frac{1}{1 + x^4}$$


```python
plot(1/(1 + x**4), (x, 0, 10))
```


```python
integrate(1/(1 + x**4), (x, 0, oo))
```




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



Of course not all integrals can be computed symbolically:


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




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



### Compute with Distributions
#### Heaviside Distribution

Heaviside function $\theta(x)$ has the following properties:
1. $\frac{d}{dx}\theta(x) = \delta(x)$
2. $\theta(x) = \begin{cases} 1 & x \geq 0\\ 0 & x < 0 \end{cases}$
3. $\frac{d}{dx}\max(x,0)=\theta(x)$


```python
Heaviside(9)
```




$\displaystyle 1$




```python
Heaviside(-9)
```




$\displaystyle 0$



There is some controversy about the behavior at $x=0$, Mathematica defines $\theta(0)=1$, but Maple uses $\theta(0)=\text{undefined}$. Different application areas may have specific conventions. For example, in control theory, it is common practice to assume $\theta(0)=0$ to match the Laplace transform of a DiracDelta distribution.


```python
Heaviside(0,)
```




$\displaystyle \theta\left(0\right)$



To specify the value of Heaviside at $x=0$, a second argument can be given. Omit this 2nd argument or pass None to recover the default behavior.


```python
Heaviside(0, S.Half)
```




$\displaystyle \frac{1}{2}$




```python
x = symbols("x ")
y = Heaviside(x)
plot(y, (x, -1, 1))
```

#### Dirac Distribution

DiracDelta function $\delta(x)$ has the following properties:

1. $\frac{d}{dx}\theta(x)=\delta(x)$
2. $\int_{a-\epsilon}^{a+\epsilon}\delta(x-a)f(x)=f(a)$
3. $\dots$


```python
dydx = diff(y)
dydx
```




$\displaystyle \delta\left(x\right)$




```python
integrate(dydx,(x,-.1,.1))
```




$\displaystyle 1$




```python
integrate(dydx,(x,.1,.2))
```




$\displaystyle 0$




```python
x = symbols("x ")
y = DiracDelta(x)
plot(y, (x, -1, 1))
```


```python
DiracDelta(0)
```




$\displaystyle \delta\left(0\right)$



## Ordinary differential equations (ODEs)

Example: $y''+4y=e^x$


```python
x = symbols("x ")
f = Function('f')
dsolve(Derivative(f(x), x, x) + 4*f(x) - exp(x), f(x))
```




$\displaystyle f{\left(x \right)} = C_{1} \sin{\left(2 x \right)} + C_{2} \cos{\left(2 x \right)} + \frac{e^{x}}{5}$



We can also add initial conditions, e.g. $y(0)=1$, $y'(0)=0$:


```python
dsolve(Derivative(f(x), x, x) + 4*f(x) - exp(x), f(x), ics={f(0): 1, f(x).diff(x).subs(x, 0): 0})
```




$\displaystyle f{\left(x \right)} = \frac{e^{x}}{5} - \frac{\sin{\left(2 x \right)}}{10} + \frac{4 \cos{\left(2 x \right)}}{5}$



**Exercise**: Solve the ODE of [Bernoulli type](https://en.wikipedia.org/wiki/Bernoulli_differential_equation) $y^3 - x^2 + x y^2 y' = 0$ with initial condition $y(1)=1$.


```python
x = Symbol('x')
f = Function('f')
dsolve(f(x)**4-x**2+x*f(x)**2*Derivative(f(x), x), f(x), ics={f(1): 1})
```




$\displaystyle f{\left(x \right)} = 1 + \left(x - 1\right)^{2} - \frac{5 \left(x - 1\right)^{3}}{3} + \frac{23 \left(x - 1\right)^{4}}{12} - \frac{26 \left(x - 1\right)^{5}}{15} + O\left(x^{6}\right)$



### Systems of ODEs

Example: $f'=3f + 4g$, $g'=-4f+3g$


```python
x = Symbol('x')
f = Function('f')
g = Function('g')
eq = (Eq(Derivative(f(x),x),3*f(x)+4*g(x)),Eq(Derivative(g(x),x),-4*f(x)+3*g(x)))
dsolve(eq)
```




    [Eq(f(x), C1*exp(3*x)*sin(4*x) + C2*exp(3*x)*cos(4*x)),
     Eq(g(x), C1*exp(3*x)*cos(4*x) - C2*exp(3*x)*sin(4*x))]



Again we can add initial conditions, e.g. $f(0)=0$, $g(0)=1$:


```python
S = dsolve(eq, ics={f(0): 0, g(0): 1})
```


```python
S[0].rhs
```




$\displaystyle e^{3 x} \sin{\left(4 x \right)}$




```python
S[1].rhs
```




$\displaystyle e^{3 x} \cos{\left(4 x \right)}$



## Code generation

SymPy allows one to [generate directly compilable and executable code](https://docs.sympy.org/latest/modules/codegen.html) in a variety of different programming languages from SymPy expressions. In addition, there are functions that generate Python importable objects that can evaluate SymPy expressions very efficiently.

- [tutorial: https://www.sympy.org/scipy-2017-codegen-tutorial/](https://www.sympy.org/scipy-2017-codegen-tutorial/)

### Code Printer

Code printer convert SymPy expressions into the target language. A list of supported languages is [here](https://docs.sympy.org/latest/modules/printing.html). 


```python
pretty_print(S[1].rhs)
```

     3⋅x         
    ℯ   ⋅cos(4⋅x)



```python
pprint(S[1].rhs)
```

     3⋅x         
    ℯ   ⋅cos(4⋅x)



```python
print_latex(S[0])
```

    f{\left(x \right)} = e^{3 x} \sin{\left(4 x \right)}



```python
latex(S[0])
```




    'f{\\left(x \\right)} = e^{3 x} \\sin{\\left(4 x \\right)}'




```python
print_ccode(S[1].rhs)
```

    exp(3*x)*cos(4*x)



```python
ccode(S[1].rhs)
```




    'exp(3*x)*cos(4*x)'




```python
mathematica_code(S[1].rhs)
```




    'Exp[3*x]*Cos[4*x]'




```python
jscode(S[1].rhs)
```




    'Math.exp(3*x)*Math.cos(4*x)'




```python
fcode(S[1].rhs)
```




    '      exp(3.0d0*x)*cos(4.0d0*x)'




```python
julia_code(S[1].rhs)
```




    'exp(3*x).*cos(4*x)'




```python
octave_code(S[1].rhs)
```




    'exp(3*x).*cos(4*x)'




```python
rust_code(S[1].rhs)
```




    '(3*x).exp()*(4*x).cos()'




```python
pycode(S[1].rhs)
```




    'math.exp(3*x)*math.cos(4*x)'




```python
print_python(S[1].rhs)
```

    x = Symbol('x')
    e = exp(3*x)*cos(4*x)


**Exercise**: Plotting SymPy functions with JavaScript

In this call we load the [Chart.js library](https://www.chartjs.org):

%%javascript
# # require.config({
# #     paths: {
# #         'chartjs': '//cdnjs.cloudflare.com/ajax/libs/Chart.js/2.6.0/Chart'
# #     }
# # });

Here is some JavaScript to do the plotting. The code will take two mathematical expressions in JavaScript and plot the functions:


```python
js_template = """\
require(['chartjs'], function(chartjs){{
function f(x) {{
    return {top_function}
}};
function g(x) {{
    return {bottom_function}
}};
function linspace(a,b,n) {{
    // From: https://gist.github.com/joates/6584908
    if(typeof n === "undefined") n = Math.max(Math.round(b-a)+1,1);
    if(n<2) {{ return n===1?[a]:[]; }}
    var i,ret = Array(n);
    n--;
    for(i=n;i>=0;i--) {{ ret[i] = (i*b+(n-i)*a)/n; }}
    return ret;
}}
var ctx = document.getElementById("{chart_id}");
var data = {{
    labels: linspace(-7.5, 7.5, 500),
    datasets: [{{
        label: "top",
        function: f,
        borderColor: "rgba(75, 192, 192, 1)",
        data: [],
        fill: false,
        lineTension: 0,
    }},
    {{
        label: "bottom",
        function: g,
        borderColor: "rgba(153, 102, 255, 1)",
        data: [],
        fill: false,
        lineTension: 0,
    }}]
}};
chartjs.Chart.pluginService.register({{
    beforeInit: function(chart) {{
        var data = chart.config.data;
        for (var i = 0; i < data.datasets.length; i++) {{
            for (var j = 0; j < data.labels.length; j++) {{
                var fct = data.datasets[i].function,
                    x = data.labels[j],
                    y = fct(x);
                data.datasets[i].data.push(y);
            }}
        }}
    }}
}});
var myBarChart = new chartjs.Chart(ctx, {{
    type: 'line',
    data: data,
    options: {{
        scales: {{
            yAxes: [{{
                ticks: {{
                    beginAtZero:true
                }}
            }}]
        }}
    }}
}});
}});
element.append("<canvas id='{chart_id}' width='400'></canvas>");\
"""
```


```python
print(js_template.format(top_function='*** TODO ***', 
                         bottom_function='*** TODO ***', 
                         chart_id='*** TODO ***'))
```

    require(['chartjs'], function(chartjs){
    function f(x) {
        return *** TODO ***
    };
    function g(x) {
        return *** TODO ***
    };
    function linspace(a,b,n) {
        // From: https://gist.github.com/joates/6584908
        if(typeof n === "undefined") n = Math.max(Math.round(b-a)+1,1);
        if(n<2) { return n===1?[a]:[]; }
        var i,ret = Array(n);
        n--;
        for(i=n;i>=0;i--) { ret[i] = (i*b+(n-i)*a)/n; }
        return ret;
    }
    var ctx = document.getElementById("*** TODO ***");
    var data = {
        labels: linspace(-7.5, 7.5, 500),
        datasets: [{
            label: "top",
            function: f,
            borderColor: "rgba(75, 192, 192, 1)",
            data: [],
            fill: false,
            lineTension: 0,
        },
        {
            label: "bottom",
            function: g,
            borderColor: "rgba(153, 102, 255, 1)",
            data: [],
            fill: false,
            lineTension: 0,
        }]
    };
    chartjs.Chart.pluginService.register({
        beforeInit: function(chart) {
            var data = chart.config.data;
            for (var i = 0; i < data.datasets.length; i++) {
                for (var j = 0; j < data.labels.length; j++) {
                    var fct = data.datasets[i].function,
                        x = data.labels[j],
                        y = fct(x);
                    data.datasets[i].data.push(y);
                }
            }
        }
    });
    var myBarChart = new chartjs.Chart(ctx, {
        type: 'line',
        data: data,
        options: {
            scales: {
                yAxes: [{
                    ticks: {
                        beginAtZero:true
                    }
                }]
            }
        }
    });
    });
    element.append("<canvas id='*** TODO ***' width='400'></canvas>");


SymPy functions now can be plotted by filling in the two missing expressions and the id in the above boilerplate code. Then call the JavaScript display function on that code: 


```python
from IPython.display import Javascript
```


```python
x=symbols('x')
f1=sin(x)
f2=cos(x)
Javascript(js_template.format(top_function=jscode(f1), 
                              bottom_function=jscode(f2), 
                              chart_id='sincos'))
```




    <IPython.core.display.Javascript object>



### Lamdify SymPy expressions

Here we convert SymPy expression evaluation to Python Lambda functions. SymPy evaluation is arbitrary precision and therefore very slow. In many cases it may be usefull to first create a Lambda function and then evaluate the function numerically.

Lets create a SymPy expression an record the time for the symbolic evaluation:


```python
y, y = symbols('x y')
expr = 3*x**2 + log(x**2 + y**2 + 1)
expr
```




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




```python
expr.subs({x: 17, y: 42}).evalf()
```




$\displaystyle 874.627544390489$




```python
%timeit expr.subs({x: 17, y: 42}).evalf()
```

    296 µs ± 17.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)


Now we create a python lambda function and compare the time for numeric evaluation.


```python
import math
f = lambda x, y: 3*x**2 + math.log(x**2 + y**2 + 1)
f(17,42)
```




    874.6275443904885




```python
%timeit f(17,42)
```

    2.07 µs ± 158 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)


Unsuprisingly, the numeric evaluation is much faster. However, the advantave is deminished by the time we spend coding the python lambda function. Here comes in `lambdify`, that generates python code and execute `eval` on it.  


```python
%timeit lambdify([x,y], expr, modules=['math'])(17,42)
```

    3.1 ms ± 364 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)



```python
g=lambdify([x,y], expr, modules=['math'])
g(17,42)
```




    874.6275443904885




```python
%timeit g(17,42)
```

    2.59 µs ± 317 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)


The code generation uses numpy by default (without `modules` parameter). 


```python
g=lambdify([x,y], expr)
g(np.linspace(17,18,5),42)
```




    array([874.62754439, 900.31920442, 926.38590757, 952.82765322,
           979.64444076])



We can get a specific function signature from `lambdify`. In the following example, the function expects a tripple as third argument.


```python
z = z1, z2, z3 = symbols('z:3')
expr = x*y*(z1+z2+z3)
expr
```




$\displaystyle x y \left(z_{0} + z_{1} + z_{2}\right)$




```python
lambdify([x,y,z],expr)(1,2,[3,4,5])
```




    24



**Exercise**: Plot $z(x,y)=\frac{df(x,y)}{dxdy}$ from $f(x,y)=3 x^{2} + \log{\left(x^{2} + y^{2} + 1 \right)}$ 


```python
y, y = symbols('x y')
expr = 3*x**2 + log(x**2 + y**2 + 1)
expr
```




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




```python
d2fdxdy = expr.diff(x,y)
xplot = np.outer(np.linspace(-5,5,100),np.ones(100))
yplot = xplot.T
zplot = lambdify([x, y],d2fdxdy)(xplot,yplot)
```


```python
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
%matplotlib notebook
```


```python
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot_surface(xplot,yplot,zplot,cmap=plt.cm.coolwarm)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('$%s$'% latex(d2fdxdy))
```


    <IPython.core.display.Javascript object>









    Text(0.5, 0, '$- \\frac{4 x y}{\\left(x^{2} + y^{2} + 1\\right)^{2}}$')




```python
%matplotlib notebook
```

The [**competitive Lotka–Volterra equations**](https://en.wikipedia.org/wiki/Competitive_Lotka–Volterra_equations) are a simple model of the population dynamics of species competing for some common resource. They are a system of differential equations that describe the dynamics of a pair of species that interact, one as a predator and the other as a prey.
The coupled first-order non-linear differential equations are given by

$$\frac{dN_1}{dt}=\frac{N_{1} r_{1} \left(K_{1} - N_{1} - N_{2} \alpha\right)}{K_{1}}$$

$$\frac{dN_2}{dt}=\frac{N_{2} r_{2} \left(K_{2} - N_{1} \beta - N_{2}\right)}{K_{2}}$$

where
- $N_i$ = population size of species $i$
- $t$ = time
- $r_i$ = species $i$ growth rate
- $K_i$ = environment carrying capacity
- $\alpha,\beta$ = inter-species interaction parameters.

If species 2 is the predator and species 1 is the prey, then $\alpha$ must be positive (since the predator eats the prey reducing its population) and $\beta$ must be negative since the predator needs the prey to get food.

Define the variables:


```python
# We split it into two definitions to remind us that
# a set of variables is for species 1 and the other one for species 2

# Species 1
r1, K1, N1, alpha = symbols('r_1 K_1 N_1 alpha')

# Species 2
r2, K2, N2, beta = symbols('r_2 K_2 N_2 beta')
```

We can now define the differential equations in sympy


```python
dN1_dt = r1 * N1 * (K1 - N1 - alpha * N2) / K1
dN1_dt
```




$\displaystyle \frac{N_{1} r_{1} \left(K_{1} - N_{1} - N_{2} \alpha\right)}{K_{1}}$




```python
dN2_dt = r2 * N2 * (K2 - N2 - beta * N1) / K2
dN2_dt
```




$\displaystyle \frac{N_{2} r_{2} \left(K_{2} - N_{1} \beta - N_{2}\right)}{K_{2}}$



Usually one of the first things we do when analyzing a dynamical system is to find the steady state, i.e. finding when the derivative with respect to time equals to zero.

We define the equality using the `Eq`:


```python
ss_N1 = Eq(dN1_dt, 0)
ss_N1
```




$\displaystyle \frac{N_{1} r_{1} \left(K_{1} - N_{1} - N_{2} \alpha\right)}{K_{1}} = 0$




```python
ss_N2 = Eq(dN2_dt, 0)
ss_N2
```




$\displaystyle \frac{N_{2} r_{2} \left(K_{2} - N_{1} \beta - N_{2}\right)}{K_{2}} = 0$



Solve for the species population:


```python
ss_N1_sol = solve(ss_N1, N1)
ss_N2_sol = solve(ss_N2, N2)

print('The steady state solution for both species are:')
ss_N1_sol, ss_N2_sol
```

    The steady state solution for both species are:





    ([0, K_1 - N_2*alpha], [0, K_2 - N_1*beta])



So we find that there are two steady state solutions:
1. The trivial solutions in which both species are extinct $N_1=N_2=0$
2. A solution in which both species maintain their non-zero populations.

We can use sympy's solve function to solve simultaneously for both species. The solution is returned as a dictionary where each key is the name of the respective variable.


```python
ss_dict = solve((Eq(N1, ss_N1_sol[1]), Eq(N2, ss_N2_sol[1])),(N1, N2))
ss_dict
```




    {N_1: (-K_1 + K_2*alpha)/(alpha*beta - 1),
     N_2: (K_1*beta - K_2)/(alpha*beta - 1)}



We can use scipy's numerical solver to run a simulation for this system of equations.

**Note**: we use `lambdify` to numerically evaluate the diff. equations


```python
# Generate 'lambdify' functions to numerically evaluate the diff. equations
dN1_dt_num = lambdify([N1, N2, r1, K1, alpha], dN1_dt)
dN2_dt_num = lambdify([N2, N1, r2, K2, beta], dN2_dt)


def lotka_voltera_rhs(n1_n2, t, *p):
    '''
    Set up the right hand side (rhs) function for the system 
    (necessary step to feed sympy's odeint function).
    '''
    # unpack the variables
    n1, n2 = n1_n2
    
    # unpack the parameters
    r_1, k1, a, r_2, k2, b = p
    
    # Define derivatices rom lambdified functions
    n1dot = dN1_dt_num(n1, n2, r_1, k1, a)
    n2dot = dN2_dt_num(n2, n1, r_2, k2, b)

    return np.array([n1dot, n2dot])
```

Now let's define the parameters and numerically solve the system


```python
from scipy.integrate import odeint
# let's set the parameter values.
r_1, r_2 = [0.5, 0.2]
k1, k2 = [100, 50]
a, b = [1.2, 0.2]
args = (r_1, k1, a, r_2, k2, b)

# Initial population size
n1_n2_0 = [2, 2]

# Time points we want to consider
t = np.linspace(0, 100, 1000)

# Now let's use odeint to solve the function
n1_n2 = odeint(lotka_voltera_rhs, n1_n2_0, t, args=args)
```

We will now plot the solution


```python
for i, N in enumerate([N1, N2]):
    plt.plot(t, n1_n2[:, i], label=r'${0:s}$'.format(latex(N)))
    plt.axhline(y=ss_dict[N].subs({K1: k1, K2: k2, alpha: a, beta: b}),
                linestyle='--',
                label=r'${0:s}$ steady state'.format(latex(N)))

plt.legend(loc=0, fontsize=15)
plt.xlabel('time (units consistent with $r_1$ and $r_2$)')
plt.ylabel('Population size')
```


    <IPython.core.display.Javascript object>









    Text(0, 0.5, 'Population size')



The [**generalised Lotka-Volterra equations**](https://en.wikipedia.org/wiki/Generalized_Lotka–Volterra_equation) model the dynamics of the populations $x_1, x_2, \dots, x_n$ of $n$ biological species ([lecture](https://stefanoallesina.github.io/Sao_Paulo_School/intro.html#basic-formulation)). Together, these populations can be considered as a vector $\mathbf{x}$. They are a set of ordinary differential equations given by

$$\frac{dx_i}{dt}=x_i f_i(\mathbf{x})$$

where the vector $\mathbf{f}$  is given by

$$\mathbf{f}=\mathbf{r}+\mathbf{A}\mathbf{x}$$

where $\mathbf{r}$ is a vector and $\mathbf{A}$ is a matrix known as the [community matrix](https://en.wikipedia.org/wiki/Community_matrix).


```python
y0 = Function('y0')
y1 = Function('y1')
y2 = Function('y2')
t = symbols('t')
C1, C2, C3 = symbols('C1, C2, C3')
b0,b1,b2 = symbols('b0:3')

y = Matrix([y0(t), y1(t), y2(t)])
B = Matrix([b0,b1,b2])

ode = (y.diff(t) - y).T.tolist()[0]
soln = dsolve(ode, y.T.tolist()[0])
soln
```




    [Eq(y0(t), C1*exp(t)), Eq(y1(t), C2*exp(t)), Eq(y2(t), C3*exp(t))]




```python
p1=plot(soln[0].subs({C1:1}).rhs,(t,0,10),show=False)
p2=plot(soln[1].subs({C2:2}).rhs,(t,0,10),show=False)
p3=plot(soln[2].subs({C3:3}).rhs,(t,0,10),show=False)
p1.append(p2[0])
p1.append(p3[0])
p1.show()
```


    <IPython.core.display.Javascript object>






Can this be written more general for any $n$?

79ddb23ecae2a72795185e79b44ef463d94835c5
ipynb
Jupyter Notebook
notebooks/SymbolicMathematics.ipynb
OleBo/MathSo
1f9fa0492d467c0bb0479768c503eee7723ae777
MIT
2022-01-08T23:52:57.000Z
2022-01-08T23:52:57.000Z
notebooks/SymbolicMathematics.ipynb
OleBo/MathSo
1f9fa0492d467c0bb0479768c503eee7723ae777
MIT
notebooks/SymbolicMathematics.ipynb
OleBo/MathSo
1f9fa0492d467c0bb0479768c503eee7723ae777
MIT
2022-01-02T21:14:30.000Z
2022-01-02T21:14:30.000Z
Qwen/Qwen-72B
1. YES
2. YES