# Diferential Equation Solver and Grapher on Sympy
### Instructions
- Change dot_x and dot_y for any function f(x,y) and g(x,y) you like
- For plotting, set the plot boundaries 
- For the trajectory, set the starting point and time 


```python
# Bulk Imports and setup
from IPython.display import display, Math, Latex
import matplotlib.pyplot as plt
import numpy as np

import sympy as sm
from sympy.utilities.lambdify import lambdify, implemented_function
x, y, lamda = sm.symbols('x, y, lamda')

def move_sympyplot_to_axes(p, ax):
    backend = p.backend(p)
    backend.ax = ax
    backend._process_series(backend.parent._series, ax, backend.parent)
    backend.ax.spines['right'].set_color('none')
    backend.ax.spines['bottom'].set_position('zero')
    backend.ax.spines['top'].set_color('none')
    plt.close(backend.fig)
```


```python
# Defining the equations via sympy
dot_x = x-2*y +x**2 -y**2            # f(x,y)
dot_y = -2*x +y +3*x**2-3*y**2           # g(x,y)
display(Math('\dot{x} = '+ sm.latex(dot_x)))
display(Math('\dot{y} = '+ sm.latex(dot_y)))
```


$\displaystyle \dot{x} = x^{2} + x - y^{2} - 2 y$



$\displaystyle \dot{y} = 3 x^{2} - 2 x - 3 y^{2} + y$



```python
### A matrix from dot(X)=AX (if system is linear)
# We contruct the vector <f(x,y), g(x,y)>.
A_vector = sm.Matrix([dot_x, dot_y])
# Jacobian for the A_vector. In a linear system this gives the A matrix from dot(X)=AX
A_matrix = A_vector.jacobian(([x, y]))
display(Math('\mathbf{A} = '+ sm.latex(A_matrix)))
```


$\displaystyle \mathbf{A} = \left[\begin{matrix}2 x + 1 & - 2 y - 2\\6 x - 2 & 1 - 6 y\end{matrix}\right]$



```python
# Setting equations to zero and computing Fixed Points
dot_x_Equal, dot_y_Equal = sm.Eq(dot_x, 0), sm.Eq(dot_y, 0)
equilibria = sm.solve((dot_x_Equal, dot_y_Equal), x, y, dict=True)
#equilibria = [{x: 0, y: 0}]
```


```python
# Gives the value of the linear approximation for the Second Lyapunov Method
# And the characteristic Polinomial
for eq in equilibria:
    if sm.im(eq[x])==0 and sm.im(eq[y]) ==0: # Excludes imaginary equilibria
        matrix_subbed = A_matrix.subs(eq)
        display(Math("A\Big|_{(" +sm.latex(eq[x]) + "," + sm.latex(eq[y])+ ")} = "+
                     sm.latex(matrix_subbed)))
        display(Math("p(\lambda)\Big|_{(" +sm.latex(eq[x]) + "," + sm.latex(eq[y])+ ")}="+
                     sm.latex(matrix_subbed.charpoly()))) 
```


$\displaystyle A\Big|_{(0,0)} = \left[\begin{matrix}1 & -2\\-2 & 1\end{matrix}\right]$



$\displaystyle p(\lambda)\Big|_{(0,0)}=\operatorname{PurePoly}{\left( \lambda^{2} - 2 \lambda - 3, \lambda, domain=\mathbb{Z} \right)}$



$\displaystyle A\Big|_{(\frac{7}{8},\frac{5}{8})} = \left[\begin{matrix}\frac{11}{4} & - \frac{13}{4}\\\frac{13}{4} & - \frac{11}{4}\end{matrix}\right]$



$\displaystyle p(\lambda)\Big|_{(\frac{7}{8},\frac{5}{8})}=\operatorname{PurePoly}{\left( \lambda^{2} + 3, \lambda, domain=\mathbb{Z} \right)}$



```python
# Calculates eigenvals and eigenvects for each eq
# The list is of the form [{eq}, [eigenthing1, eigenthing1]]
# Each eigenthing is [(eigenval, multiplicity, [eigenvec]), …]
eigen_things = [[eq, A_matrix.subs(eq).eigenvects()]
              for eq in equilibria]
# Showing equilibrium and eigenthings
for _, eq in enumerate(equilibria):
    if sm.im(eq[x])==0 and sm.im(eq[y]) ==0: # Excludes imaginary equilibria
        display(Math("\mathrm{Equilibrium}: (" + sm.latex(eq[x])+","+sm.latex(eq[y])+")"))
        # eigen_things[#eq][eigenthings][#pair][0=value, 1=mult, 2=eigenvect]
        for pair in eigen_things[_][1]:
                display(Math("\mathrm{Eigenvalue}:"+ sm.latex(pair[0]) +
                             "\quad \mathrm{Multiplicity}=" + sm.latex(pair[1]) + "\quad\mathrm{Eigenvector}: "+ sm.latex(pair[2][0])))

```


$\displaystyle \mathrm{Equilibrium}: (0,0)$



$\displaystyle \mathrm{Eigenvalue}:-1\quad \mathrm{Multiplicity}=1\quad\mathrm{Eigenvector}: \left[\begin{matrix}1\\1\end{matrix}\right]$



$\displaystyle \mathrm{Eigenvalue}:3\quad \mathrm{Multiplicity}=1\quad\mathrm{Eigenvector}: \left[\begin{matrix}-1\\1\end{matrix}\right]$



$\displaystyle \mathrm{Equilibrium}: (\frac{7}{8},\frac{5}{8})$



$\displaystyle \mathrm{Eigenvalue}:- \sqrt{3} i\quad \mathrm{Multiplicity}=1\quad\mathrm{Eigenvector}: \left[\begin{matrix}\frac{11}{13} - \frac{4 \sqrt{3} i}{13}\\1\end{matrix}\right]$



$\displaystyle \mathrm{Eigenvalue}:\sqrt{3} i\quad \mathrm{Multiplicity}=1\quad\mathrm{Eigenvector}: \left[\begin{matrix}\frac{11}{13} + \frac{4 \sqrt{3} i}{13}\\1\end{matrix}\right]$



```python
# Ploting  parameters
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1,1,1)
# ax.set_title("Sistema diferencial")

# Setting Plot Boundaries (needed of sympy)
x_min, y_min, x_max, y_max,  = -.2, -.2, 1.2, 1.1
# For plotting the trajectory
time, dt = 60, 0.1
initial_pairs = [(.85, .55), (.5, .5)] #List of the form [(x_0, y_0)]
# Choose if normalization of vectors is applied
normalize = True

#################################
# Plot fixed points
for point in equilibria:
    # Excludes complex equilibria
    if sm.im(point[x])==0 and sm.im(point[y]) == 0:
        ax.plot(point[x],point[y],"red", marker = "o", markersize = 10.0)
    
# To plot we translate from sympy to numpy
x_movement = lambdify((x, y), dot_x) 
y_movement = lambdify((x, y), dot_y)

# define a grid and compute direction at each point
x_line = np.linspace(x_min, x_max, 20)
y_line = np.linspace(y_min, y_max, 20)

# Plot initial point and trayectory
for pair in initial_pairs: 
    x_0 = pair[0]
    y_0 = pair[1]
    ax.plot(x_0, y_0, "black", marker = "*", markersize = 7.0)
    # Plot a singular trayectory
    x_list = [x_0]
    y_list = [y_0]
    for i in range(time):
        x_list.append(x_list[i] + (x_movement(x_list[i], y_list[i])) * dt)
        y_list.append(y_list[i] + (y_movement(x_list[i], y_list[i])) * dt)
    ax.plot(x_list, y_list, color="black")


# Quiverplot
# create a grid and set the dir vecs
X, Y = np.meshgrid(x_line, y_line)               
x_dir, y_dir = x_movement(X,Y), y_movement(X,Y)  
# prepare the colour and normalization factor
norm = np.sqrt(x_dir**2 + y_dir**2)

#Normalize
if normalize:
    x_dir /= norm
    y_dir /= norm

plt.quiver(X, Y, 
           x_dir, 
           y_dir,
           norm,
           pivot='mid')

# Plot Nullclines
p1 = sm.plot_implicit(dot_x_Equal, (x, x_min, x_max), (y, y_min, y_max), show= False)
p2 = sm.plot_implicit(dot_y_Equal, (x, x_min, x_max), (y, y_min, y_max), show= False)

move_sympyplot_to_axes(p1, ax)
move_sympyplot_to_axes(p2, ax)

plt.show()

```

Note: if "Result too large error", just decrese the time 

# Credit Where Credit is due
- http://systems-sciences.uni-graz.at/etextbook/sw2/phpl_python.html
- https://stackoverflow.com/questions/60325325/putting-together-plots-of-matplotlib-and-sympy
- Dynamical Systems with Applications Using Python - Stephen Lynch (2018)

eb8b186407080b7a1a1364406fd09b8cde4bc21d
ipynb
Jupyter Notebook
diffeq.ipynb
ndharari/diffeq
0130b774fac40c5c1865db4fbba201c3ea042f56
MIT
diffeq.ipynb
ndharari/diffeq
0130b774fac40c5c1865db4fbba201c3ea042f56
MIT
diffeq.ipynb
ndharari/diffeq
0130b774fac40c5c1865db4fbba201c3ea042f56
MIT
Qwen/Qwen-72B
1. YES
2. YES