# # ODE-Analysis.py: A basic class + main() to study and visualize the properties of linear and non-linear ODEs # Used in 15-382, Collective Intelligence # Author: Gianni A. Di Caro, CMU-Q, 2018 # Free to use acknowledging the source # # # Overview: # The class considers the case of a generic two-dimensional autonomous system of ODEs # The system can be either linear or non-linear. # # A linear dynamical model is fully specified through a model matrix A (2x2). # A non-linear model requires the definition of two model-specific equations. # A linear model can be always solved analytically and completely studied in terms of its stability. # A non-linear model can present difficulties finding general solutions in analytic form, and it # can be necessary to use a numeric integration of system's equations. # Moreover, in order to study the stability of the critical points, the study of the linear system # obtained from the linearization of the non-linear system near the critical points is normally # useful, if not necessary. # # Based on these basic observations, the class provides a number of methods to define, solve, and # study the properties of both linear and non-linear two-dimensional autonomous systems. # The code is based on the use of SymPy, an extensive library for symbolic mathematics, ScyPy, # a large library of mathematical tools, and Matplotlib, a scientific plotting library. # An example of use is provide in the main() function. ## Docs and tutorials on Sympy: # http://docs.sympy.org/latest/index.html # http://docs.sympy.org/latest/tutorial/index.html # ## A tutorial tailored for a course on dynamical systems: # http://www.cfm.brown.edu/people/dobrush/am33/SymPy/index.html from __future__ import division import numpy as np import matplotlib.pyplot as plt # ONLY FOR my F.. MAC OS # it doesn't find the package by itself! import sys, os #print '\n'.join(sys.path) # to check the import paths if sys.platform == 'darwin': sys.path.append(os.path.join(os.path.dirname(__file__), "/Library/Python/2.7/site-packages/sympy-1.1.1-py2.7.egg/")) from sympy import * from sympy.plotting import plot from sympy.plotting import plot_parametric from sympy.abc import * import sympy as sm from scipy.integrate import odeint import sys class DynamicalSystem(): def __init__(self): print( "---------- Declaration of general symbols and parameters ----------------") # In the case of linear system, the general form of the dynamical model is a 2x2 matrix A: # dx1/dt = a11*x1 + a12*x2 = f1(x1,x2) # dx2/dt = a21*x1 + a22*x2 = f1(x1,x2) # Solutions are in the form x1(t), x2(t) # (f1, f2) is the vector field, where the functions f1 and f2 are linear # # For a non-linear model: # dx1/dt = f1(x1, x2) # dx2/dt = f2(x1, x2) # Solutions are in the form x1(t), x2(t) # (f1, f2) is the vector field, where the functions f1 and f2 are generically non-linear # Declaration of the symbols that are used to derive a solution and manipulate the dynamical model self.x1, self.x2, self.x3 = sm.symbols("x1 x2 x3", cls = Function, Function = True) # solution functions x1(t), x2(t) self._x1, self._x2, self._x3 = sm.symbols('self._x1, self._x2, self._x3', negative=False) # state variables (x1, x2) self.t = sm.symbols('t') # time parameter self.C1, self.C2 = sm.symbols('C1 C2') # C1 and C2 are integration constants # Linear model matrix A: an empty 2x2 matrix of type sm.Matrix is created # self.A = zeros(2) # A number of general parameters need to be specified, that are used # at plotting time, to show vector fields and vector flows in the phase portraits. # Each dynamical model is better shown using model-specific # parameter values for the ranges of the (x1,x2) points, as well as of the time. # The values used below are just general reasonable default values that can / should # be changed when the model is instantiated in linear_model() and nonlinear_model() methods. # # The time_span parameter is used to plot trajectories from t=0 to t=timespan. # The selected ranges for x1 and x2 (xmin / xmax) are used to define for which # (x1, x2) points the behavior of the system is shown. # The parameters numx1 and numx2 are the number of plotted points along each axis. # self.time_span = 1 self.trajectory_pts = 200 self.x1min = -1.5; self.x1max = 1.5; self.x1step = 0.5; self.x2min = -1.0; self.x2max = 1.0; self.x2step = 0.5; self.numx1 = 20; self.numx2 = 20; init_printing() # useful if some advanced printing properties of SymPy are going to be used ... print( "") def linear_model(self, model_type="spiral"): # print( "---------- Linear dynamical model definition ----------------") # # Calling this method allows to set up a specific linear dynamical system model, both # in terms of its numeric coefficient matrix A, and as a symbolic system of differential equations. # A list of a few example models are available and selectable using argument "model_type". # If a new model is needed, it must be added below, in the body of the "else" statement. # For each model instance, are specified: the coefficients of the model matrix A, # the time/space parameters for visualization of the phase portraits (feel free to change them) # if model_type == "spiral": # Spiral #a11 = -1; a12 = -10; a21 = 10; a22 = -1 self.A = sm.Matrix([[-1, -10], [10, -1]]) self.time_span = 5 self.trajectory_pts = 250 self.x1min = -1.5; self.x1max = 1.5; self.x1step = 0.5; self.x2min = -1.0; self.x2max = 1.0; self.x2step = 0.5; self.numx1 = 20; self.numx2 = 20; elif model_type == "center": self.A = sm.Matrix([[0, 1], [-1, -0]]) self.time_span = 10 self.trajectory_pts = 200 self.x1min = -1.5; self.x1max = 1.5; self.x1step = 0.25; self.x2min = -1.0; self.x2max = 1.0; self.x2step = 0.25; self.numx1 = 20; self.numx2 = 20; elif model_type == "saddle": self.A = sm.Matrix([[-4, -3], [2, 3]]) self.time_span = 0.5 self.trajectory_pts = 300 self.x1min = -5; self.x1max = 5; self.x1step = 1; self.x2min = -5; self.x2max = 5; self.x2step = 1; self.numx1 = 20; self.numx2 = 20; elif model_type == "node": # self.A = sm.Matrix([[3, 1], [0, 3]]) # self.time_span = 1 # self.trajectory_pts = 300 # self.x1min = -5; self.x1max = 5; self.x1step = 1; # self.x2min = -5; self.x2max = 5; self.x2step = 1; # self.numx1 = 20; self.numx2 = 20; self.A = sm.Matrix([[-2, 0], [1, -4]]) self.time_span = 0.5 self.trajectory_pts = 300 self.x1min = -5; self.x1max = 5; self.x1step = 1; self.x2min = -5; self.x2max = 5; self.x2step = 1; self.numx1 = 20; self.numx2 = 20; elif model_type == "star": self.A = sm.Matrix([[-1, 0], [0, -1]]) self.time_span = 2 self.trajectory_pts = 400 self.x1min = -5; self.x1max = 5; self.x1step = 1; self.x2min = -5; self.x2max = 5; self.x2step = 1; self.numx1 = 30; self.numx2 = 30; elif model_type == "improper node": self.A = sm.Matrix([[-1, 2], [0, -1]]) self.time_span = 2 self.trajectory_pts = 400 self.x1min = -5; self.x1max = 5; self.x1step = 1; self.x2min = -5; self.x2max = 5; self.x2step = 1; self.numx1 = 30; self.numx2 = 30; else: #model_type == "new" # here write your specific model and study it! print( "Add the definition of the new dynamical model in the model() method!") sys.exit(-2) # Using the above model specification (the matrix A) a symbolic representation of # the system of differential equations is created. # # diff(x1(t), t) represents dx1/dt, # diff(x2(t), t) represents dx2/dt, # Together, they fully specifies the vector field associated to A # self.f1 = sm.Eq(diff(self.x1(t), self.t), self.A[0]*self.x1(self.t) + self.A[1]*self.x2(self.t)) self.f2 = sm.Eq(diff(self.x2(t), self.t), self.A[2]*self.x1(self.t) + self.A[3]*self.x2(self.t)) # # if we would remove the 'self' parts a more readable form would be obtained: #f1 = sm.Eq(diff(x1(t), t), A[0]*x1(t) + A[1]*x2(t)) #f2 = sm.Eq(diff(x2(t), t), A[2]*x1(t) + A[3]*x2(t)) print( "ODE linear system model:\n", self.f1, "\n", self.f2) # Symbolic system of equations representing only the vector field # This is used to find the critical points # Note that in this case there's no dependence on t (that is necessary for describing the flows) # The symbolic variables are _x1, _x2, that represent the state variables # self.X1 = self.A[0] * self._x1 + self.A[1] * self._x2 self.X2 = self.A[2] * self._x1 + self.A[3] * self._x2 # The symbolic model definitions X1, X2, are "lambdifyied" # to efficiently compute the values and treat them as a numeric function. # This is not used in this method, but rather in the vector_field() method # self.X1lambdaLin = lambdify((self._x1, self._x2), self.X1, "numpy") self.X2lambdaLin = lambdify((self._x1, self._x2), self.X2, "numpy") print( "") def nonlinear_model_competing_species(self): print( "---------- Non-Linear dynamical model definition: Competing species ----------------") # # Definition of a non-linear model of a dynamical system. The model is fully specified # both as a symbolic system of equations, and as a numeric system of functions. # Since each non-linear model is specified by a different set of equations, it's not possible # to provide a general form, but only specific guiding examples. # The symbolic model could be of limited use in the case of non-linear systems since # it may be too diffult to use it for finding an analytic solution (using SymPy). # In any case, this method is also used to assign a number of parameters that are necessary # for the study of the non-linear model. # # Example model: Competing species # parameters of the model (these could be also defined as symbols, to be free parameters of solutions) # self.growth1 =1 self.growth2 = 3/4 self.saturation1 = 1 self.saturation2 = 1 self.inference21 = 1/3 #1 self.inference12 = 1/2 # Symbolic system of differential equations. The formalism is analogous to the linear case. # self.f1 = sm.Eq(diff(self.x1(t), self.t), self.x1(self.t) * (self.growth1 - self.saturation1 * self.x1(self.t) - self.inference21 * self.x2(self.t) )) self.f2 = sm.Eq(diff(self.x2(t), self.t), self.x2(self.t) * (self.growth2 - self.saturation1 * self.x2(self.t) - self.inference12 * self.x1(self.t) )) print( "[Competing populations] ODE non-linear system model:\n", self.f1, "\n", self.f2) # Symbolic system of equations representing only the vector field # This is used to derive the Jacobian and to find the critical points # Note that in this case there's no dependence on t (that is necessary for describing the flows) # self.X1 = self._x1 * (self.growth1 - self.saturation1 * self._x1 - self.inference21 * self._x2) self.X2 = self._x2 * (self.growth2 - self.saturation2 * self._x2 - self.inference12 * self._x1) # The symbolic model definitions X1, X2, are "lambdifyied" # to efficiently compute the values and treat them as a numeric function. # This is not used in this method, but rather in the vector_field() method # self.X1lambdaNonLin = lambdify((self._x1, self._x2), self.X1, "numpy") self.X2lambdaNonLin = lambdify((self._x1, self._x2), self.X2, "numpy") # parameters for studying and plotting / representing the model self.time_span = 40 self.trajectory_pts = 250 self.x1min = 0; self.x1max = 1; self.x1step = 0.05; self.x2min = 0; self.x2max = 1; self.x2step = 0.05; self.numx1 = 20; self.numx2 = 20; # # nullclines: I should do it automatically, but for the time being is manual :-) self.x1nc_x1 = np.linspace(self.x1min, self.x1max, 100) self.x1nc_x2 = (self.growth1 - self.saturation1 * self.x1nc_x1) / self.inference21 self.x2nc_x1 = np.linspace(self.x2min, self.x2max, 100) self.x2nc_x2 = (self.growth2 - self.inference12 * self.x2nc_x1) / self.saturation2 print( "") def nonlinear_model_eval(self, X, t): # The function computes the value of the model for the given point # X corresponds to a specific point in the phase space, X = (x1, x2) # The function exploit the lambdifyied version of the vector fields # defined in the model method # #print( "---------- Numeric model evaluation ----------------") # f1 = self.X1lambdaNonLin(X[0], X[1]) f2 = self.X2lambdaNonLin(X[0], X[1]) #print( "Model value in (x1=%5.2f, x2=%5.2f): (%5.2f, %5.2f)\n" % (X[0], X[1], f1, f2)) #print( "") return [f1, f2] def general_solution_symbolic(self, f1, f2): # print( "---------- Finding general solution (works for linear systems) ----------------") # # Solution of the system, ODEsol is tuple [x1(t), x2(t)] in symbolic form # self.ODEsol = sm.dsolve((f1, f2)) print( "Solutions of the ODE system:\n", self.ODEsol[0], "\n", self.ODEsol[1]) print( "") def numeric_solutions_and_system_flows(self): # print( "---------- Finding and plotting solutions in numeric form (for non linear systems) --------------") # # plt.figure("Flows (numeric solutions)") plt.xlabel('$x_1$') plt.ylabel('$x_2$') plt.xlim([self.x1min, self.x1max]) plt.ylim([self.x2min, self.x2max]) for x1_0 in np.arange(self.x1min, self.x1max, self.x1step): for x2_0 in np.arange(self.x2min, self.x2max, self.x2step): # this defines the range for the time parameter and the number of points: trajectories will # represent the range [0 : timespan] using trajectory_pts points equally spaced # timespan = np.linspace(0, self.time_span, self.trajectory_pts) Xsol = odeint(self.nonlinear_model_eval, [x1_0, x2_0], timespan) plt.plot(Xsol[:,0], Xsol[:,1], 'b-') # trajectory path, blue #plt.plot([Xsol[0,0]], [Xsol[0,1]], 'o') # trajectory start, circle #plt.plot([Xsol[-1,0]], [Xsol[-1,1]], 's') # trajectory end, square for (cx, cy) in self.critical_pts: #print( "[%5.2f, %5.2f]" % (cx, cy)) plt.plot(cx, cy, marker='o', markersize=7, color="black") plt.savefig('flows-numeric.png') plt.ioff() plt.gcf().show() # keep the window open print( "") def eigenspaces(self, M): # The argument M is a matrix representing a linear system, likely self.A # print( "---------- Eigenanalysis ----------------") # # Stability study requires eigenvalues and eigenvectors # print( "Determinant: %8.3f" % M.det()) #print( "Null space", M.nullspace()) print( "Eigenvalues:", M.eigenvals()) # eigenvalues can be real or complex number, print them in a different format accordingly for (e,m) in M.eigenvals().items(): if sm.im(e) == 0: print( "\teigenvalue: %5.2f multiplicity: %d" % (e,m)) else: if sm.im(e) < 0: print( "\teigenvalue: %5.2f %5.2fi multiplicity: %d" % (sm.re(e), sm.im(e), m)) else: print( "\teigenvalue: %5.2f + %5.2fi multiplicity: %d" % (sm.re(e), sm.im(e), m)) print( "Eigenvectors:") for (e,m,v) in M.eigenvects(): if sm.im(e) == 0: print( "\t[%5.2f,%5.2f] <- eigenvalue %5.2f" % (v[0][0], v[0][1], e)) else: print( "\t[(%5.2f,%5.2f), (%5.2f,%5.2f)] <- eigenvalue %5.2f + %5.2fi" % (sm.re(v[0][0]), sm.im(v[0][0]), sm.re(v[0][1]), sm.im(v[0][1]), sm.re(e), sm.im(e))) print( "") def critical_points(self, X1, X2, system_type): # The function arguments, X1 and X2 are the vector fields in symbolic form, # that must have _x1, _x2 as variables # print( "---------- Critical points ----------------") # # This method finds the critical points. In the case of a linear system a linear-specific # method could be effectively used. In the non-linear case a more general method is needed. # setting equations to zero X1Eq = sm.Eq(X1, 0) X2Eq = sm.Eq(X2, 0) # compute critical points if system_type == "linear": self.critical_pts = sm.linsolve( (X1Eq, X2Eq), self._x1, self._x2 ) else: self.critical_pts = sm.solve( (X1Eq, X2Eq), self._x1, self._x2 ) for xy in self.critical_pts: if len(xy) == 2: print( "[%5.2f, %5.2f]" % (xy[0], xy[1])) else: print( xy ) return self.critical_pts def def_Jacobian_matrix(self, VF1, VF2): # The arguments are the two components of a (non-linear) vector field # print( "---------- Definition of Jacobian matrix ----------------") # F = sm.Matrix([VF1, VF2]) #print( "F: ", F) self.Jsymb = F.jacobian([self._x1, self._x2]) #print( "Jacobian matrix:\n", self.J) print( "Jacobian matrix:\n", "\t", self.Jsymb[0, 0], "\t\t", self.Jsymb[0, 1], "\n\t", self.Jsymb[1, 0], "\t\t", self.Jsymb[1, 1]) print( "") # def Jacobian_eval(self, J, [x1coord, x2coord]): def Jacobian_eval(self, J, x1x2coord): # value of the Jacobian in the specified point x1x2coord = (x1coord, x2coord) # print( "---------- Jacobian evaluation ----------------") j_x1x2 = J.subs([(self._x1, x1x2coord[0]), (self._x2, x1x2coord[1])]) print( "Jacobian at coordinates (%5.2f, %5.2f)" % (x1x2coord[0], x1x2coord[1])) print( "\t %5.2f \t %5.2f\n\t %5.2f \t %5.2f" % (j_x1x2[0, 0], j_x1x2[0, 1], j_x1x2[1, 0], j_x1x2[1, 1])) print( "") return j_x1x2 def vector_field(self, VectorFieldX1, VectorFieldX2, nullclines = True): # The arguments VectorFieldX1 and VectorFieldX2 are lambdifyed versions of the vector field # of the model. They are set up in the critical_points() method, but are passed as an argument # leaving more flexibility defining them # # print( "---------- Phase portrait: vector field ----------------") # # To generate the phase portrait, we need to use the derivatives dx1/dt and # dx2/dt at t=0 on a grid over the range of values for x1 and x2 we are # interested in. That is, we need to compute the vector field F=(u,v): # a vector (u,v) at each (x1, x2) that shows the tangent direction of the # (unique) trajectory passing by the point. _x1points = np.linspace(self.x1min, self.x1max, self.numx1) _x2points = np.linspace(self.x2min, self.x2max, self.numx2) # a regular meshgrid is defined over the selected ranges in the phase space X1grid, X2grid = np.meshgrid(_x1points, _x2points) # init the vector field arrays u and v to correct dimensions and to zeros # (un, vn) are used to represent a normalized vector field (all vectors have same length) # u, v = np.zeros(X1grid.shape), np.zeros(X2grid.shape) un, vn = np.zeros(X1grid.shape), np.zeros(X2grid.shape) # compute the value of the vector field over the selected grid points using the lambdifyied # versions of the model equations # for i in range(self.numx1): for j in range(self.numx2): x1coord = X1grid[i, j] # from grid to cartesian coordinates x2coord = X2grid[i, j] # unnormalized (regular) vector field: vectors have length proportional to the norm u[i,j] = VectorFieldX1(x1coord, x2coord) v[i,j] = VectorFieldX2(x1coord, x2coord) # normalized: all vectors have same length (maybe it's more readable, but less info) norm = np.sqrt(u[i,j]**2.0 + v[i,j]**2.0) if norm > 0: un[i,j] = u[i,j] / np.sqrt(u[i,j]**2.0 + v[i,j]**2.0) vn[i,j] = v[i,j] / np.sqrt(u[i,j]**2.0 + v[i,j]**2.0) else: un[i,j] = 0.0 vn[i,j] = 0.0 # plot the arrows representing the vector field (u,v) (and the normalized version) in # the selected grid points X1grid, X2grid using the function quiver() # plt.figure("Vector field") plt.xlabel('$x_1$') plt.ylabel('$x_2$') plt.xlim([self.x1min, self.x1max]) plt.ylim([self.x2min, self.x2max]) plt.quiver(X1grid, X2grid, u, v, color='r') for (cx, cy) in self.critical_pts: #print( "[%5.2f, %5.2f]" % (cx, cy)) plt.plot(cx, cy, marker='o', markersize=7, color="black") plt.savefig('vector-field.png') plt.ioff() plt.gcf().show() # keep the window open #plt.ion() #plt.show() plt.figure("Vector field (uniform vector length)") plt.xlabel('$x_1$') plt.ylabel('$x_2$') plt.xlim([self.x1min, self.x1max]) plt.ylim([self.x2min, self.x2max]) plt.quiver(X1grid, X2grid, un, vn, color='b') for (cx, cy) in self.critical_pts: #print( "[%5.2f, %5.2f]" % (cx, cy)) plt.plot(cx, cy, marker='o', markersize=7, color="black") plt.savefig('vector-field-uniform.png') # nullclines if nullclines == True: plt.plot(self.x1nc_x1, self.x1nc_x2) plt.plot(self.x2nc_x1, self.x2nc_x2) plt.savefig('nullclines.png') plt.ioff() plt.gcf().show() print( "") def vector_flow_linear(self, X1sol, X2sol): # print( "---------- Phase portrait: vector flow ----------------") # # Selection of a set of specific trajectories by assigning initial conditions # In turn, this allows to set the values for the integration constants in the general sol # The multiple parametric trajectories are plotted in the same phase portrait # flows_fig = plt.figure("Flows") ax = flows_fig.add_subplot(111) # creates one plot that will include all selected trajectories for x1_0 in np.arange(self.x1min, self.x1max, 2*self.x1step): for x2_0 in np.arange(self.x2min, self.x2max, 2*self.x2step): # trajectory_constants = sm.solve( (self.ODEsol[0].subs(self.t,0).subs(self.x1(0), x1_0), # self.ODEsol[1].subs(self.t,0).subs(self.x2(0), x2_0)), # {self.C1, self.C2}) trajectory_constants = sm.solve( (X1sol.subs(self.t,0).subs(self.x1(0), x1_0), X2sol.subs(self.t,0).subs(self.x2(0), x2_0)), {self.C1, self.C2}) # # subs method substitutes symbol t with 0, and x1(0), x2(0) with the selected # initial conditions, the result is used to solve with respect to the integration # constants C1 and C2 (previously declared as symbols) show_parametric_equations = False if show_parametric_equations: print( "Integration constants from initial conditions (%5.2f, %5.2f):" % (x1_0, x2_0)) print( trajectory_constants) print( "") # Solution (in symbolic form) for the selected trajectory # two (parametric) equations dependending on t # x1trajectory = sm.expand(self.ODEsol[0].rhs.subs(trajectory_constants)) x2trajectory = sm.expand(self.ODEsol[1].rhs.subs(trajectory_constants)) if show_parametric_equations: print( "Solution trajectory in the (x1, x2) phase space:") print( "x1(t)=", x1trajectory, " x2(t)=", x2trajectory) print( "") # found trjectories are in symbolic form, in order to use them to generate a set of points # to be plotted, they need to be lambdifyied, in order to have numeric type of functions x1t_numeric = lambdify(self.t, x1trajectory, "numpy") x2t_numeric = lambdify(self.t, x2trajectory, "numpy") # this defines the range for the time parameter and the number of points: trajectories will # represent the range [0 : timespan] using trajectory_pts points equally spaced timespan = np.linspace(0, self.time_span, self.trajectory_pts) x1t_pts = x1t_numeric(timespan) x2t_pts = x2t_numeric(timespan) ax.plot(x1t_pts, x2t_pts) plt.xlabel('$x_1$') plt.ylabel('$x_2$') plt.ioff() plt.gcf().show() # keep the window open flows_fig.savefig('./vector-flow.png') print( "") #import prey_predator if __name__ == '__main__': linear = False #True nonlinear = True #True # # Linear system analysis # if linear == True: DSLin = DynamicalSystem() DSLin.linear_model("star") DSLin.general_solution_symbolic(DSLin.f1, DSLin.f2) DSLin.eigenspaces(DSLin.A) DSLin.critical_points(DSLin.X1, DSLin.X2, "linear") DSLin.vector_field(DSLin.X1lambdaLin, DSLin.X2lambdaLin, False) DSLin.vector_flow_linear(DSLin.ODEsol[0], DSLin.ODEsol[1]) input("Press Enter to continue ...") plt.close('all') del DSLin # Non-Linear system analysis # if nonlinear == True: DSNonLin = DynamicalSystem() DSNonLin.nonlinear_model_competing_species() #DSNonLin.general_solution_symbolic(DSNonLin.f1, DSNonLin.f2) critical_points = DSNonLin.critical_points(DSNonLin.X1, DSNonLin.X2, "non-linear") DSNonLin.vector_field(DSNonLin.X1lambdaNonLin, DSNonLin.X2lambdaNonLin, False) DSNonLin.numeric_solutions_and_system_flows() DSNonLin.def_Jacobian_matrix(DSNonLin.X1, DSNonLin.X2) for cp in critical_points: J = DSNonLin.Jacobian_eval(DSNonLin.Jsymb, [cp[0], cp[1]]) DSNonLin.eigenspaces(J) input("Press Enter to exit...") plt.close('all')