# # ODE-Analysis.py: A basic class + main() to study and visualize properties of linear and non-linear ODEs # Used in 15-382, Collective Intelligence, Spring 2018 # Author: Gianni A. Di Caro, CMU-Q, Jan 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. # # Usage: # ----- # An example of use is provide in the __main()__ function. Two separate models, one linear, and one # non-linear (Lotka-Volterra) are considered. Just run the python code. # ** In order to study a new linear model: add the model (the matrix A), in the method linear_model() # # ** In order to study a new non-linear model XYZ: add a new method, possibly named nonlinear_model_XYZ(), # following the way the lotka-volterra model is specified in nonlinear_model_prey_predator() # # In both cases, in the same methods, adjust the visualization parameters for # the vector fields and trajectory flows # ## 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 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 = sm.symbols("x1 x2", cls = Function, Function = True) # solution functions x1(t), x2(t) self._x1, self._x2 = sm.symbols('self._x1, self._x2', 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 = 300 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 = 250 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 = 10 self.trajectory_pts = 400 self.x1min = -5; self.x1max = 5; self.x1step = 0.75; self.x2min = -5; self.x2max = 5; self.x2step = 0.75; self.numx1 = 30; self.numx2 = 30; 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; 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. Once defined the matrix A above, # the rest of the method doesn't need any further edit. # # 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_prey_predator(self): # # 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. # print "---------- Non-Linear dynamical model definition: Prey-Predator ----------------" # # Prey-Predator, Lotka-Volterra model # parameters of the model # self.growth_prey = 1 self.death_pred = 3/4 self.interact_predprey = 1/2 self.interact_preypred = 1/4 # 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.growth_prey - self.interact_predprey * self.x2(self.t)) ) self.f2 = sm.Eq(diff(self.x2(t), self.t), self.x2(self.t) * (-self.death_pred + self.interact_preypred * 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.growth_prey - self.interact_predprey * self._x2) self.X2 = self._x2 * (-self.death_pred + self.interact_preypred * 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 = 10 self.trajectory_pts = 120 self.x1min = 0.0; self.x1max = 9; self.x1step = 0.5; self.x2min = 0.0; self.x2max = 7; self.x2step = 0.5; self.numx1 = 20; self.numx2 = 20; # 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 exploits 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, 2*self.x1step): for x2_0 in np.arange(self.x2min, self.x2max, 2*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 plt.ioff() plt.gcf().show() # keep the window open plt.savefig('./flows-non-linear.png') 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. # Setting equations to zero X1Eq = sm.Eq(X1, 0) X2Eq = sm.Eq(X2, 0) # Compute critical points if system_type == "linear": critical_pts = sm.linsolve( (X1Eq, X2Eq), self._x1, self._x2 ) else: critical_pts = sm.solve( (X1Eq, X2Eq), self._x1, self._x2 ) for (cx, cy) in critical_pts: print "[%5.2f, %5.2f]" % (cx, cy) print "" return 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)): # value of the Jacobian in the specified point (x1coord, x2coord) # print "---------- Jacobian evaluation ----------------" j_x1x2 = J.subs([(self._x1, x1coord), (self._x2, x2coord)]) print "Jacobian at coordinates (%5.2f, %5.2f)" % (x1coord, x2coord) 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): # 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') plt.savefig('vector-field.png') plt.ioff() plt.gcf().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') plt.savefig('vector-field-uniform.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 plt.xlim([self.x1min, self.x1max]) plt.ylim([self.x2min, self.x2max]) 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-linear.png') print "" if __name__ == '__main__': linear = True #False # nonlinear = True #False # # Linear system analysis # if linear == True: DSLin = DynamicalSystem() DSLin.linear_model("saddle") # set up the model DSLin.general_solution_symbolic(DSLin.f1, DSLin.f2) # find solutions x(t) DSLin.critical_points(DSLin.X1, DSLin.X2, "linear") # non necessary for non degenerate A DSLin.eigenspaces(DSLin.A) # eigenanalysis of the origin DSLin.vector_field(DSLin.X1lambdaLin, DSLin.X2lambdaLin) # compute and plot vector field DSLin.vector_flow_linear(DSLin.ODEsol[0], DSLin.ODEsol[1]) # plot trajectory flows raw_input("Press Enter to continue ...") plt.close('all') del DSLin # Non-Linear system analysis # if nonlinear == True: DSNonLin = DynamicalSystem() DSNonLin.nonlinear_model_prey_predator() # set up the model #DSNonLin.nonlinear_model_competing_species() #DSn.general_solution_symbolic(DS.f1, DS.f2) # attempt finding general sol critical_points = DSNonLin.critical_points(DSNonLin.X1, DSNonLin.X2, "non-linear") # find critical points DSNonLin.vector_field(DSNonLin.X1lambdaNonLin, DSNonLin.X2lambdaNonLin) # compute and plot vector field DSNonLin.numeric_solutions_and_system_flows() # numeric solutions and trajectory flows DSNonLin.def_Jacobian_matrix(DSNonLin.X1, DSNonLin.X2) # define general Jacobian matrix for cp in critical_points: J = DSNonLin.Jacobian_eval(DSNonLin.Jsymb, [cp[0], cp[1]]) # evaluate Jacobian in each critical point DSNonLin.eigenspaces(J) # eigenanalysis of the Jacobian raw_input("Press Enter to exit...") plt.close('all')