CoCalc Shared FilesPhase portrait.sagews
Author: Johann Thiel
Views : 44
##################################################
# Phase portrait generator
##################################################
# 3/22 - Added functionality: The function
# phasePortrait() will now plot multiple solution
# curves when given a list of points. See the
# example below.

##################################################
# User inputs example
##################################################
# User-defined dynamical system
#f1(x1,x2) = x1-x2
#f2(x1,x2) = x1+x1*x2

# User-defined initial points for phase portrait
#points = [[1,1.2],[2,3],[-1,-1]]

# User-defined graphing window
#xRange = [-5,5]
#yRange = [-5,5]
##################################################

##################################################
# Script
##################################################
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def phasePortrait(pts,xR,yR,g1,g2):

myList = []

for point in pts:

def g(x,t):
X1 = x[0]
X2 = x[1]
G1 = g1(X1,X2)
G2 = g2(X1,X2)
return [G1,G2]

# initial conditions
x0 = point
t  = np.linspace(0, 100, 5000)

# numerical solution
soln = odeint(g, x0, t)
X1 = soln[:, 0]
X2 = soln[:, 1]
myList.append([(a,b) for a,b in zip(X1,X2) if xR[0]<= a <= xR[1] and yR[0]<= b <= yR[1]])

P = Graphics()

for L in myList:
P+= list_plot(L)

P += plot_vector_field((g1,g2),(x1,xR[0],xR[1]),(x2,yR[0],yR[1]))

show(P)

#phasePortrait(points,xRange,yRange,f1,f2)