- Email: [email protected]
Introduction to Phase Plane Analysis
C H A P T E R
14 Introduction to Phase Plane Analysis 14.1 GOAL OF THIS CHAPTER The goal of this chapter is to examine the cone and horizontal cell system using a qualitative visualization technique called phase plane analysis; this system will be discussed further in Chapter 28, “Models of the Retina.” The techniques presented here will be used again in Chapter 15, “Exploring the Fitzhugh-Nagumo Model.”
14.2 BACKGROUND In this chapter you will be studying a retinal feedback model; this model is described further in Chapter 28, “Models of the Retina.” The system is represented as follows: dC~ 1 ~ 5 ð2 C~ 2 kHÞ dt τC
ð14:1Þ
dH~ 1 ~ 5 ð2 H~ 1 CÞ dt τH
ð14:2Þ
Typical values for these parameters are τ C 5 0.025 sec, τ H 5 0.08 sec, and k 5 4. Now assume that the light intensity is L 5 10 (i.e., daylight). For your initial conditions, choose that C(0) 5 H(0) 5 0. Finally, be aware that: L C~ 5 C 2 k11
and
L H~ 5 H 2 k11
ð14:3Þ
For further details of the basic biology of this system, see Chapter 28, “Models of the Retina.” In that chapter, we will examine in more detail the model of retinal feedback between cone cells and horizontal cells of the retina, shown in Equations 14.1 and 14.2. Although the explicit solutions determined in that chapter are more informative, many more complicated systems (such as the Fitzhugh-Nagumo system presented in Chapter 15, “Exploring the Fitzhugh-Nagumo Model”) can only be qualitatively described. When we describe a system
MATLAB® for Neuroscientists. DOI: http://dx.doi.org/10.1016/B978-0-12-383836-0.00014-X
253
© 2014 Elsevier Inc. All rights reserved.
254
14. INTRODUCTION TO PHASE PLANE ANALYSIS
qualitatively, we look for steady-state values of the solutions (often called fixed points) and try to classify the dynamics of the solution that led to these steady-state values. In Chapter 28, “Models of the Retina,” we will consider the following system in more detail: dx 5x1y dt dy 5 4x 1 y dt which has eigenvalues 21 and 3 and has the solution:
ð14:4Þ
xðtÞ 5 C1 e3t 1 C2 e2t
ð14:6Þ
yðtÞ 5 2C1 e3t 2 2C2 e2t
ð14:7Þ
ð14:5Þ
This solution can be described qualitatively. If you wait long enough, then this system will approach one of two states. If C1 5 0, then: limt-N xðtÞ 5 limt-N yðtÞ 5 0
ð14:8Þ
Therefore, one says that (x, y) 5 (0, 0) is a steady-state or fixed point of the system. For C1 ¼ 6 0, then: limt-N xðtÞ 5 limt-N yðtÞ 5 N
ð14:9Þ
Therefore, the only finite steady-state solution to this system is (x, y) 5 (0, 0). Regardless of how you choose C1 and C2, there are no other steady-state values for this system. Since the initial conditions determine C1 and C2, then those initial conditions that lead to C1 5 0 will have solutions that steadily tend toward the fixed point (0, 0), while all others will steadily tend toward infinity (i.e., away from the fixed point at the origin). A fixed point with this property—that is, with some initial conditions leading to the fixed point and others leading away from it—is called a saddle point. This simple qualitative description of identifying the steady state(s) of the solution, the dynamics of what initial conditions lead to the steady state(s), and how it is reached steadily or in an oscillatory fashion can all be determined from a phase plane analysis of the system. The first step in phase plane analysis is to set up a phase plane. The axes for the plane represent the state variables characterizing the system. In the preceding example, the phase plane is constructed with y as the ordinate and x as the abscissa. Next, the -x- and y-nullclines are plotted. The x-nullcline is the curve in the x-y plane, where: dx 50 dt A similar definition applies for the y-nullcline. Intersections of these nullclines represent points where: dx dy 5 50 dt dt so x and y are no longer changing with time. In other words, these intersections represent steady-state values or fixed points of the system. Next, a vector field is constructed by assigning the following vector to every point on the x-y plane:
III. DATA ANALYSIS WITH MATLAB
255
14.2 BACKGROUND
"
dx dt
dy dt
#T
Notice that this vector field can be determined without knowing the solution to the system. Since the slope of these vectors is: dy dx dy 5 ð14:10Þ m5 dt dt dx by the chain rule, the vector field must be tangent to any solution (x, y) of the system. This allows you to use the vector field to calculate the solution of the system for any initial condition (xo,yo). Such a solution when plotted on the phase plane is called a trajectory. The phase plane, nullclines, vector field, and several trajectories are shown in Figure 14.1 for the system in Equations 14.4 and 14.5. x’ = x + y y’ = 4 x + y 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10
−8
−6
−4
Cursor position: (1.79, –15.8)
−2
0 x
2
4
The backward orbit from (–2.7, –4.8) left the computation window. Ready. The forward orbit from (1.9, –5.4) left the computation window. The backward orbit from (1.9, –5.4) left the computation window. Ready.
FIGURE 14.1 Phase plane of a linear system showing saddle node stability.
III. DATA ANALYSIS WITH MATLAB
6
8
10
256
14. INTRODUCTION TO PHASE PLANE ANALYSIS
In the figure, the nullclines are plotted as dashed lines. Notice that these nullclines intersect at the point (x, y) 5 (0, 0) indicating that this is the steady-state of the system in agreement with what was predicted by considering the explicit solutions (Equations 14.6 and 14.7). Any linear system of ordinary differential equations described by a matrix with real eigenvalues of opposite sign (recall that the eigenvalues for this system are 21 and 3) will have a saddle point at the intersection of its nullclines. If the matrix describing the linear system has real eigenvalues that are both negative, then the fixed point is called a nodal sink. The classic phase portrait of a nodal sink is shown in Figure 14.2. If the matrix describing the linear system has real eigenvalues that are both positive, then the fixed point is called a nodal source. The classic phase portrait of a nodal source is shown in Figure 14.3. Notice the difference in the direction of the arrows in the vector field in this figure. x’ = −x y’ = −3y 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10
−8
−6
Cursor position: (−9.75, 12)
−4
−2
0 x
2
4
6
The forward orbit from (−2.9, −3.8) a possible eq. pt. near (−5.6e-014, −5.1e-020). The backward orbit from (−2.9, −3.8) left the computation window. The forward orbit from (−6.1, −4.1) a possible eq. pt. near (−8.4e-014, −2e-020). The backward orbit from (−6.1, −4.1) left the computation window. Ready.
FIGURE 14.2
Phase plane of a linear system showing nodal sink stability.
III. DATA ANALYSIS WITH MATLAB
8
10
257
14.2 BACKGROUND
x’ = −x y’ = 3y 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10 −8 −6 −4 Cursor position: (−9.58, 12.1)
−2
0 x
2
4
6
8
10
The backward orbit from (−2.5, −5.1) a possible eq. pt. near (−6e-014, 8.1e-019). Ready. The forward orbit from (−4.3, −4.9) left the computation window. The backward orbit from (−4.3, −4.9) a possible eq. pt. near (−7e-014, −3.7e-020). Ready.
FIGURE 14.3 Phase plane of a linear system showing nodal source stability.
If the matrix describing the linear system has imaginary eigenvalues that have negative real parts, then the fixed point is called a spiral sink. The classic phase portrait of a spiral sink is shown in Figure 14.4. If the matrix describing the linear system has imaginary eigenvalues that have positive real parts, then the fixed point is called a spiral source. The classic phase portrait of a spiral source is shown in Figure 14.5. These five types of equilibria are collectively known as the generic equilibria. There are also five nongeneric equilibria. The most important nongeneric equilibrium is called a center. It occurs when the eigenvalues of the matrix are purely imaginary. The classic phase portrait of a center is shown in Figure 14.6.
III. DATA ANALYSIS WITH MATLAB
258
14. INTRODUCTION TO PHASE PLANE ANALYSIS
x’ = 4x − 3y y’ = 15x − 8y 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10
−8
−6
−4
Cursor position: (−11.3, 12.2)
−2
0 x
2
4
6
8
10
The backward orbit from (−0.82, −3.2) left the computation window. Ready. The forward orbit from (4.2, −4.6) a possible eq. pt. near (1.5e-014, −1e-014). The backward orbit from (4.2, −4.6) left the computation window. Ready.
FIGURE 14.4
Phase plane of a linear system showing spiral sink stability.
14.3 EXERCISES The phase portraits in the preceding section were drawn using a downloadable M-file called pplane7.m. The phase plane consists of three basic features: the nullclines intersecting at the fixed point of the system, the vector field showing how the solutions change over time, and trajectories showing how the solution approaches its steady-state from a given initial condition. The first exercise of this chapter will involve writing a simple version of pplane7. Several functions built into MATLAB will aid in coding each of the basic features of the phase plane mentioned previously. Plotting the nullclines of the system requires no more than the basic plotting commands used throughout previous chapters. Plotting the vector fields can be greatly aided by the functions meshgrid() and quiver(). The function meshgrid takes two vector arguments x and y, and returns two square matrices X and Y such that each row of X is a copy of the
III. DATA ANALYSIS WITH MATLAB
259
14.3 EXERCISES
x’ = 2x − y y’ = 2x 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10 −8 −6 −4 Cursor position: (−9.28, 12.2)
−2
0 x
2
4
6
8
10
The backward orbit from (−1.5, 6.1) a possible eq. pt. near (−6.6e-016, −1.7e-014). Ready. The forward orbit from (4.7, 4.5) left the computation window. The backward orbit from (4.7, 4.5) a possible eq. pt. near (−6.9e-015, −7.1e-015). Ready.
FIGURE 14.5 Phase plane of a linear system showing spiral source stability.
vector x, and each column of Y is a copy of the vector y. This function is useful for evaluating functions of two variables. For example, suppose you wanted to evaluate the function f(x,y) 5 x 1 y. You could do this using for loops; for example, you could type .. .. .. .. .. .. ..
x 5 [0:0.1:10]; y 5 x; for ii 5 1:length(x) for jj 5 1:length(y) f(ii,jj) 5 x(ii) 1 y(jj); end; end;
which produces the same results as the commands
III. DATA ANALYSIS WITH MATLAB
260
14. INTRODUCTION TO PHASE PLANE ANALYSIS
x’ = − y y’ = x 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10 −8 −6 −4 Cursor position: (−11.1, 12.2)
−2
0 x
2
4
6
8
10
The backward orbit from (3, −4.7) a nearly closed orbit. Ready. The forward orbit from (4.1, −8.6) a nearly closed orbit. The backward orbit from (4.1, −8.6) a nearly closed orbit. Ready.
FIGURE 14.6
.. .. .. ..
Phase plane of a linear system showing center stability
x 5 210:10; y 5 x; [X, Y] 5 meshgrid(x,y); f 5 X 1 Y;
Evaluating functions of two variables is important in this chapter because the model you wish to study (Equations 14.1 and 14.2) expresses the derivatives as functions of the two ~ By comparing the matrix f as defined in the preceding code to Equation variables: C~ and H. dx for several values of x and y. You 14.1, you see that f holds the values of the derivative dt could define a matrix g that holds the y-derivative by using the following command: .. g 5 4X 1 Y;
III. DATA ANALYSIS WITH MATLAB
261
14.3 EXERCISES
10 8 6 4 2 0 –2 –4 –6 –8 –10 –10
–8
–6
–4
–2
0
2
4
6
8
10
FIGURE 14.7 Vector field created by the quiver() command.
Once you have evaluated these derivatives using meshgrid, we can plot a vector field using the quiver command. If you have the matrices X, Y, f, and g defined as shown here, then type this command: .. quiver(X,Y,f,g); You should get the result shown in Figure 14.7. The function quiver works by plotting a vector on the plane at points (x, y) with components (f, g). You can plot the trajectory of a system given an initial condition in several ways. One method is to use a numerical solver such as the ode_euler() or RK4() functions you will write in Chapter 25, “Voltage-Gated Ion Channels,” to solve for x and y given some initial condition and then plot x versus y. Another method would be to calculate the derivatives of x and y at the initial condition, move the system a short distance in the direction indicated by the derivative, and then repeat over many time steps. Either method will work, and both can be done with no more than the basic functions introduced in Chapter 2, “MATLAB Tutorial.” EXERCISE 14.1 Write a function phase_plane(A, init) that takes a matrix and performs a phase
plane analysis for the linear system, u0 5 Au. The function should plot a phase
III. DATA ANALYSIS WITH MATLAB
262
14. INTRODUCTION TO PHASE PLANE ANALYSIS
plane with axes x and y, plot the nullclines, create the vector field, and plot the phase plane trajectory that passes through the initial condition. Finally, the program should output the type of equilibrium point, saddle point, spiral sink, etc. If the equilibrium point is nongeneric, then the program can just output nongeneric as the class type. This is quite a complicated program, so you may want to write several smaller functions
that can be called within phase_plane—for example, a separate function that will simply classify the fixed point and then another that will create a vector field, etc. Hint: The Boolean function isreal() will return 0 if the argument is not a real number and 1 if it is. This function might be useful for deciding whether or not the eigenvalues of A are real, so that the fixed point of the system can be classified.
14.4 PROJECT Use your phase_plane program to analyze the retinal model described at the beginning of this chapter. The matrix describing this linear system is: 2 3 21 2 k 6 τC τC 7 6 7 6 1 21 7 4 5 τH τH Identify what kind of behavior the fixed point exhibits. Repeat using the parameters for dim light: τ C 5 0:1 sec; τ H 5 0:5 sec; and k 5 0:5 What is the behavior of the fixed point now?
MATLAB FUNCTIONS, COMMANDS, AND OPERATORS COVERED IN THIS CHAPTER isreal quiver eig
III. DATA ANALYSIS WITH MATLAB
14 Introduction to Phase Plane Analysis 14.1 GOAL OF THIS CHAPTER The goal of this chapter is to examine the cone and horizontal cell system using a qualitative visualization technique called phase plane analysis; this system will be discussed further in Chapter 28, “Models of the Retina.” The techniques presented here will be used again in Chapter 15, “Exploring the Fitzhugh-Nagumo Model.”
14.2 BACKGROUND In this chapter you will be studying a retinal feedback model; this model is described further in Chapter 28, “Models of the Retina.” The system is represented as follows: dC~ 1 ~ 5 ð2 C~ 2 kHÞ dt τC
ð14:1Þ
dH~ 1 ~ 5 ð2 H~ 1 CÞ dt τH
ð14:2Þ
Typical values for these parameters are τ C 5 0.025 sec, τ H 5 0.08 sec, and k 5 4. Now assume that the light intensity is L 5 10 (i.e., daylight). For your initial conditions, choose that C(0) 5 H(0) 5 0. Finally, be aware that: L C~ 5 C 2 k11
and
L H~ 5 H 2 k11
ð14:3Þ
For further details of the basic biology of this system, see Chapter 28, “Models of the Retina.” In that chapter, we will examine in more detail the model of retinal feedback between cone cells and horizontal cells of the retina, shown in Equations 14.1 and 14.2. Although the explicit solutions determined in that chapter are more informative, many more complicated systems (such as the Fitzhugh-Nagumo system presented in Chapter 15, “Exploring the Fitzhugh-Nagumo Model”) can only be qualitatively described. When we describe a system
MATLAB® for Neuroscientists. DOI: http://dx.doi.org/10.1016/B978-0-12-383836-0.00014-X
253
© 2014 Elsevier Inc. All rights reserved.
254
14. INTRODUCTION TO PHASE PLANE ANALYSIS
qualitatively, we look for steady-state values of the solutions (often called fixed points) and try to classify the dynamics of the solution that led to these steady-state values. In Chapter 28, “Models of the Retina,” we will consider the following system in more detail: dx 5x1y dt dy 5 4x 1 y dt which has eigenvalues 21 and 3 and has the solution:
ð14:4Þ
xðtÞ 5 C1 e3t 1 C2 e2t
ð14:6Þ
yðtÞ 5 2C1 e3t 2 2C2 e2t
ð14:7Þ
ð14:5Þ
This solution can be described qualitatively. If you wait long enough, then this system will approach one of two states. If C1 5 0, then: limt-N xðtÞ 5 limt-N yðtÞ 5 0
ð14:8Þ
Therefore, one says that (x, y) 5 (0, 0) is a steady-state or fixed point of the system. For C1 ¼ 6 0, then: limt-N xðtÞ 5 limt-N yðtÞ 5 N
ð14:9Þ
Therefore, the only finite steady-state solution to this system is (x, y) 5 (0, 0). Regardless of how you choose C1 and C2, there are no other steady-state values for this system. Since the initial conditions determine C1 and C2, then those initial conditions that lead to C1 5 0 will have solutions that steadily tend toward the fixed point (0, 0), while all others will steadily tend toward infinity (i.e., away from the fixed point at the origin). A fixed point with this property—that is, with some initial conditions leading to the fixed point and others leading away from it—is called a saddle point. This simple qualitative description of identifying the steady state(s) of the solution, the dynamics of what initial conditions lead to the steady state(s), and how it is reached steadily or in an oscillatory fashion can all be determined from a phase plane analysis of the system. The first step in phase plane analysis is to set up a phase plane. The axes for the plane represent the state variables characterizing the system. In the preceding example, the phase plane is constructed with y as the ordinate and x as the abscissa. Next, the -x- and y-nullclines are plotted. The x-nullcline is the curve in the x-y plane, where: dx 50 dt A similar definition applies for the y-nullcline. Intersections of these nullclines represent points where: dx dy 5 50 dt dt so x and y are no longer changing with time. In other words, these intersections represent steady-state values or fixed points of the system. Next, a vector field is constructed by assigning the following vector to every point on the x-y plane:
III. DATA ANALYSIS WITH MATLAB
255
14.2 BACKGROUND
"
dx dt
dy dt
#T
Notice that this vector field can be determined without knowing the solution to the system. Since the slope of these vectors is: dy dx dy 5 ð14:10Þ m5 dt dt dx by the chain rule, the vector field must be tangent to any solution (x, y) of the system. This allows you to use the vector field to calculate the solution of the system for any initial condition (xo,yo). Such a solution when plotted on the phase plane is called a trajectory. The phase plane, nullclines, vector field, and several trajectories are shown in Figure 14.1 for the system in Equations 14.4 and 14.5. x’ = x + y y’ = 4 x + y 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10
−8
−6
−4
Cursor position: (1.79, –15.8)
−2
0 x
2
4
The backward orbit from (–2.7, –4.8) left the computation window. Ready. The forward orbit from (1.9, –5.4) left the computation window. The backward orbit from (1.9, –5.4) left the computation window. Ready.
FIGURE 14.1 Phase plane of a linear system showing saddle node stability.
III. DATA ANALYSIS WITH MATLAB
6
8
10
256
14. INTRODUCTION TO PHASE PLANE ANALYSIS
In the figure, the nullclines are plotted as dashed lines. Notice that these nullclines intersect at the point (x, y) 5 (0, 0) indicating that this is the steady-state of the system in agreement with what was predicted by considering the explicit solutions (Equations 14.6 and 14.7). Any linear system of ordinary differential equations described by a matrix with real eigenvalues of opposite sign (recall that the eigenvalues for this system are 21 and 3) will have a saddle point at the intersection of its nullclines. If the matrix describing the linear system has real eigenvalues that are both negative, then the fixed point is called a nodal sink. The classic phase portrait of a nodal sink is shown in Figure 14.2. If the matrix describing the linear system has real eigenvalues that are both positive, then the fixed point is called a nodal source. The classic phase portrait of a nodal source is shown in Figure 14.3. Notice the difference in the direction of the arrows in the vector field in this figure. x’ = −x y’ = −3y 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10
−8
−6
Cursor position: (−9.75, 12)
−4
−2
0 x
2
4
6
The forward orbit from (−2.9, −3.8) a possible eq. pt. near (−5.6e-014, −5.1e-020). The backward orbit from (−2.9, −3.8) left the computation window. The forward orbit from (−6.1, −4.1) a possible eq. pt. near (−8.4e-014, −2e-020). The backward orbit from (−6.1, −4.1) left the computation window. Ready.
FIGURE 14.2
Phase plane of a linear system showing nodal sink stability.
III. DATA ANALYSIS WITH MATLAB
8
10
257
14.2 BACKGROUND
x’ = −x y’ = 3y 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10 −8 −6 −4 Cursor position: (−9.58, 12.1)
−2
0 x
2
4
6
8
10
The backward orbit from (−2.5, −5.1) a possible eq. pt. near (−6e-014, 8.1e-019). Ready. The forward orbit from (−4.3, −4.9) left the computation window. The backward orbit from (−4.3, −4.9) a possible eq. pt. near (−7e-014, −3.7e-020). Ready.
FIGURE 14.3 Phase plane of a linear system showing nodal source stability.
If the matrix describing the linear system has imaginary eigenvalues that have negative real parts, then the fixed point is called a spiral sink. The classic phase portrait of a spiral sink is shown in Figure 14.4. If the matrix describing the linear system has imaginary eigenvalues that have positive real parts, then the fixed point is called a spiral source. The classic phase portrait of a spiral source is shown in Figure 14.5. These five types of equilibria are collectively known as the generic equilibria. There are also five nongeneric equilibria. The most important nongeneric equilibrium is called a center. It occurs when the eigenvalues of the matrix are purely imaginary. The classic phase portrait of a center is shown in Figure 14.6.
III. DATA ANALYSIS WITH MATLAB
258
14. INTRODUCTION TO PHASE PLANE ANALYSIS
x’ = 4x − 3y y’ = 15x − 8y 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10
−8
−6
−4
Cursor position: (−11.3, 12.2)
−2
0 x
2
4
6
8
10
The backward orbit from (−0.82, −3.2) left the computation window. Ready. The forward orbit from (4.2, −4.6) a possible eq. pt. near (1.5e-014, −1e-014). The backward orbit from (4.2, −4.6) left the computation window. Ready.
FIGURE 14.4
Phase plane of a linear system showing spiral sink stability.
14.3 EXERCISES The phase portraits in the preceding section were drawn using a downloadable M-file called pplane7.m. The phase plane consists of three basic features: the nullclines intersecting at the fixed point of the system, the vector field showing how the solutions change over time, and trajectories showing how the solution approaches its steady-state from a given initial condition. The first exercise of this chapter will involve writing a simple version of pplane7. Several functions built into MATLAB will aid in coding each of the basic features of the phase plane mentioned previously. Plotting the nullclines of the system requires no more than the basic plotting commands used throughout previous chapters. Plotting the vector fields can be greatly aided by the functions meshgrid() and quiver(). The function meshgrid takes two vector arguments x and y, and returns two square matrices X and Y such that each row of X is a copy of the
III. DATA ANALYSIS WITH MATLAB
259
14.3 EXERCISES
x’ = 2x − y y’ = 2x 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10 −8 −6 −4 Cursor position: (−9.28, 12.2)
−2
0 x
2
4
6
8
10
The backward orbit from (−1.5, 6.1) a possible eq. pt. near (−6.6e-016, −1.7e-014). Ready. The forward orbit from (4.7, 4.5) left the computation window. The backward orbit from (4.7, 4.5) a possible eq. pt. near (−6.9e-015, −7.1e-015). Ready.
FIGURE 14.5 Phase plane of a linear system showing spiral source stability.
vector x, and each column of Y is a copy of the vector y. This function is useful for evaluating functions of two variables. For example, suppose you wanted to evaluate the function f(x,y) 5 x 1 y. You could do this using for loops; for example, you could type .. .. .. .. .. .. ..
x 5 [0:0.1:10]; y 5 x; for ii 5 1:length(x) for jj 5 1:length(y) f(ii,jj) 5 x(ii) 1 y(jj); end; end;
which produces the same results as the commands
III. DATA ANALYSIS WITH MATLAB
260
14. INTRODUCTION TO PHASE PLANE ANALYSIS
x’ = − y y’ = x 10 8 6 4
y
2 0 −2 −4 −6 −8 −10 −10 −8 −6 −4 Cursor position: (−11.1, 12.2)
−2
0 x
2
4
6
8
10
The backward orbit from (3, −4.7) a nearly closed orbit. Ready. The forward orbit from (4.1, −8.6) a nearly closed orbit. The backward orbit from (4.1, −8.6) a nearly closed orbit. Ready.
FIGURE 14.6
.. .. .. ..
Phase plane of a linear system showing center stability
x 5 210:10; y 5 x; [X, Y] 5 meshgrid(x,y); f 5 X 1 Y;
Evaluating functions of two variables is important in this chapter because the model you wish to study (Equations 14.1 and 14.2) expresses the derivatives as functions of the two ~ By comparing the matrix f as defined in the preceding code to Equation variables: C~ and H. dx for several values of x and y. You 14.1, you see that f holds the values of the derivative dt could define a matrix g that holds the y-derivative by using the following command: .. g 5 4X 1 Y;
III. DATA ANALYSIS WITH MATLAB
261
14.3 EXERCISES
10 8 6 4 2 0 –2 –4 –6 –8 –10 –10
–8
–6
–4
–2
0
2
4
6
8
10
FIGURE 14.7 Vector field created by the quiver() command.
Once you have evaluated these derivatives using meshgrid, we can plot a vector field using the quiver command. If you have the matrices X, Y, f, and g defined as shown here, then type this command: .. quiver(X,Y,f,g); You should get the result shown in Figure 14.7. The function quiver works by plotting a vector on the plane at points (x, y) with components (f, g). You can plot the trajectory of a system given an initial condition in several ways. One method is to use a numerical solver such as the ode_euler() or RK4() functions you will write in Chapter 25, “Voltage-Gated Ion Channels,” to solve for x and y given some initial condition and then plot x versus y. Another method would be to calculate the derivatives of x and y at the initial condition, move the system a short distance in the direction indicated by the derivative, and then repeat over many time steps. Either method will work, and both can be done with no more than the basic functions introduced in Chapter 2, “MATLAB Tutorial.” EXERCISE 14.1 Write a function phase_plane(A, init) that takes a matrix and performs a phase
plane analysis for the linear system, u0 5 Au. The function should plot a phase
III. DATA ANALYSIS WITH MATLAB
262
14. INTRODUCTION TO PHASE PLANE ANALYSIS
plane with axes x and y, plot the nullclines, create the vector field, and plot the phase plane trajectory that passes through the initial condition. Finally, the program should output the type of equilibrium point, saddle point, spiral sink, etc. If the equilibrium point is nongeneric, then the program can just output nongeneric as the class type. This is quite a complicated program, so you may want to write several smaller functions
that can be called within phase_plane—for example, a separate function that will simply classify the fixed point and then another that will create a vector field, etc. Hint: The Boolean function isreal() will return 0 if the argument is not a real number and 1 if it is. This function might be useful for deciding whether or not the eigenvalues of A are real, so that the fixed point of the system can be classified.
14.4 PROJECT Use your phase_plane program to analyze the retinal model described at the beginning of this chapter. The matrix describing this linear system is: 2 3 21 2 k 6 τC τC 7 6 7 6 1 21 7 4 5 τH τH Identify what kind of behavior the fixed point exhibits. Repeat using the parameters for dim light: τ C 5 0:1 sec; τ H 5 0:5 sec; and k 5 0:5 What is the behavior of the fixed point now?
MATLAB FUNCTIONS, COMMANDS, AND OPERATORS COVERED IN THIS CHAPTER isreal quiver eig
III. DATA ANALYSIS WITH MATLAB