Preface

---

This tutorial was made solely for the purpose of education and it was designed for students taking Applied Math 0330. It is primarily for students who have very little experience or have never used Mathematica and programming before and would like to learn more of the basics for this computer algebra system. As a friendly reminder, don't forget to clear variables in use and/or the kernel. The Mathematica commands in this tutorial are all written in bold black font, while Mathematica output is in normal font.

Finally, you can copy and paste all commands into your Mathematica notebook, change the parameters, and run them because the tutorial is under the terms of the GNU General Public License (GPL). You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. The tutorial accompanies the textbook Applied Differential Equations. The Primary Course by Vladimir Dobrushkin, CRC Press, 2015; http://www.crcpress.com/product/isbn/9781439851043

We need the following result that follows from the linearization procedure.

Bifurcation

---

Many of the equations that we have examined have a parameter, which means that we actually have a family of differential equations. For example, the logistic equation has two parameters, r and N. In this section we will investigate how the solutions of a differential equation vary as we change the value of a parameter or parameters.

Consider an autonomous differential equation depending on a parameter k: \( {\text d}y /{\text d}t = f(y;k) . \) We say that a bifurcation occurs at the parameter value k = k₀ if there is a change in the qualitative nature of the families of solutions when the parameter varies in a neighborhood of k₀.

A bifurcation diagram for a parameterized family of autonomous differential equations depending on a parameter k, \( {\text d}y /{\text d}t = f(y;k) , \) is a diagram in the ky-plane that summarizes the change in qualitative behavior of the equation as the parameter k is varied. The word bifurcate means “to divide into two parts or branches.” Hence, a bifurcation diagram shows us at what parameter values additional solutions emerge (or disappear). The bifurcation diagram is constructed by plotting the parameter value k against all corresponding equilibrium values \( y^{\ast} . \) Typically, k is plotted on the horizontal axis and critical points y^* on the vertical axis. A "curve" of sinks is indicated by a solid line and a curve of sources is indicated by a dashed line. The method of constructing bifurcation diagrams should become evident as we study the following examples.

We need to make an observation that bifurcation occurs at the point \( \left( k_{0} , y^{\ast} \right) \) only when

A bifurcation does not happen:

• if y^* is a sink with \( f\left( y^{\ast} \right) < 0 , \) then for all k sufficiently close to k₀ the critical point y^* remains a sink because the inequality \( f\left( y^{\ast} \right) < 0 \) holds in some neighborhood of the equilibrium point;
• if y^* is a source with \( f\left( y^{\ast} \right) > 0, \) then for all k sufficiently close to k₀ the critical point y^* remains a source because the inequality \( f\left( y^{\ast} \right) > 0 \) holds in some neighborhood of the equilibrium point.

Example: Logistic Models with Harvesting Example: autonomous equation

Example: Let us find the bifurcation values of the one-parameter family

If \( g_{k} (y) = g(y;k) = y\left( y-3\right)^2 + k , \) then \( g'_{k} (y) = \left( y-3\right)^2 + y \left( y-3\right) = 2\,y^2 -9y +9 . \) The roots of \( g'_{k} (y) = 0 \) are y =3 and y = 3/2. In order for k to be a bifurcation value, we must have \( g_{k} (3) = k = 0 \) or Therefore, the given equation has two bifurcation points k = -8/27 and k = 0. Upon solving the cube equation \( y\left( y-3\right)^2 + k =0 , \) we obtain the real critical point By plotting the above function, we get the bifurcation diagram for this one-parameter family:

 

I. Saddle-node bifurcation

---

A saddle-node bifurcation is a local bifurcation in which two (or more) critical points (or equilibria) of a differential equation (or a dynamic system) collide and annihilate each other. Saddle-node bifurcations may be associated with hysteresis and catastrophes.

Consider the slope function \( f(x, \alpha ) , \) where α is a control parameter. In this example, we use α instead of k because we need clearly distinguish a parameter (Greek letter α) from the variable (Latin letter y); we will return to k after rescaling. Then the solutions of the autonomous differential equation \( x' = f(x, \alpha ) \) and critical points (if any) will also depend on this parameter α. If we demand \( f(x, \alpha ) =0 \) and \( \partial_x f(x, \alpha ) =0, \) we get two equations in two unknowns, and in general there will be a zero-dimensional solution set consisting of points \( \left( x_c , \alpha_c \right) . \) Let’s expand \( f(x, \alpha ) \) in the vicinity of such a point \( \left( x_c , \alpha_c \right) : \)

where we keep terms of lowest order in the deviations δx and δα. If we now rescale \( u \equiv \sqrt{B/A} \left( x- x_c \right) , \quad k \equiv \left( \alpha - \alpha_c \right) , \) and \( \tau = Ax , \) we obtain, neglecting the higher order terms, the ‘normal form’ or generic model of the saddle-node bifurcation:

The right-hand side of this differential equation is the function \( f(u;k) = k+ u^2 \) that depends continuously on the parameter k. If we graph this function as a function of u for different values of k to determine equilibria, we see that for k < 0, there are two equilibria, which disappear when k becomes positive.

We can also graph the equilibria as functions of k. This figure is another example of a bifurcation diagram. For k < 0, there are two fixed points -- one stable (\( u^{\ast} = - \sqrt{-k} \) ) and one unstable (\( u^{\ast} = + \sqrt{-k} \) ). At k = 0, these two nodes coalesce and annihilate each other. (The point u^* = 0 is half-stable precisely at k = 0.) For k > 0, there are no longer any fixed points in the vicinity of u = 0. In the left panel of Figure, we show the flow in the extended (k,u) plane. The unstable and stable nodes annihilate at k = 0.

It is possible to calculate the time it takes a solution starting at \( -\infty \) at τ=0 to reach \( \infty \) at τ = T. We can find T by solving the differential equation \( \dot{u} = k + u^2 \) with the conditions \( u(0) = -\infty \) and \( u(T) = \infty . \) We separate variables and compute the following definite integrals:

The integral on the right-hand side is equal to T. The integral on the left-hand side needs more work. We will use the substitution \( x= u/\sqrt{k} \) and hence \( \sqrt{k}\,{\text d}x = {\text d}u . \) Upon factoring k in the denominator, we get Therefore, We thus find that the time it takes the solution to go from \( -\infty \) at τ=0 to reach \( \infty \) at τ = T is \( \pi / \sqrt{k} . \) This shows that as k approaches 0 from the right, it takes the solution longer and longer. The slowdown occurs when the solution passes through u = 0.

 

II. Transcritical bifurcation

---

A transcritical bifurcation is a particular kind of local bifurcation when stability of critical points changes as the parameter is varied. In other words, the unstable fixed point becomes stable and vice versa. A critical point under a transcritical bifurcation is never destroyed, it just interchanges its stability with another critical point.

Another situation which arises frequently is the transcritical bifurcation. Consider the equation \( \dot{y} = f(y) \) in the vicinity of a fixed point y^*.

We rescale \( u \equiv \beta \left( y - y^{\ast} \right) \) with \( \beta = - \frac{1}{2}\, f'' \left( y^{\ast} \right) \) and define \( k = f' \left( y^{\ast} \right) \) as the control parameter, to obtain, to order u² Note that the sign of the u² term can be reversed relative to the others by sending \( u \,\mapsto -u \) and \( k \,\mapsto -k . \) The right-hand side of this differential equation is the function \( f(u;k) = k\, u - u^2 \) that depends continuously on the parameter k. If we graph this function as a function of u for different values of k to determine equilibria, we see that there are two equilibria. These equilibria change stability as k changes from a negative to a positive value.

Example. Consider a crude model of a laser threshold. Let n be the number of photons in the laser cavity, and N the number of excited atoms in the cavity. The dynamics of the laser are approximated by the equations

Here G is the gain coefficient and k the photon decay rate. N₀ is the pump strength, and α is a numerical factor. The first equation tells us that the number of photons in the cavity grows with a rate GN-k; gain is proportional to the number of excited atoms, and the loss rate is a constant cavity-dependent quantity (typically through the ends, which are semi-transparent). The second equation says that the number of excited atoms is equal to the pump strength minus a term proportional to the number of photons (since the presence of a photon means an excited atom has decayed). Putting them together, which exhibits a transcritical bifurcation at pump strength N₀ = k/G. For N₀ < k/G the system acts as a lamp; for N₀ > k/G the system acts as a laser.

 

III. Pitchfork bifurcation

---

The pitchfork bifurcation is commonly encountered in systems in which there is an overall parity symmetry (\( u\,\mapsto -u . \) ). There are two classes of pitchfork: supercritical and subcritical. The normal form or generic model of the supercritical bifurcation is

which has fixed points at u^* = 0 and \( u^{\ast} = \pm \sqrt{k} . \) Thus, the situation is as depicted in Figure. For k < 0, there is a single stable fixed point at u^* = 0. For k > 0, u^* = 0 is unstable, and flanked by two stable fixed points at \( u^{\ast} = \pm \sqrt{k} . \)

If we send \( u\,\mapsto -u , \quad k\,\mapsto -k , \quad \mbox{and} \quad t\,\mapsto -t , \) we obtain the subcritical pitchfork, depicted in Figure. The normal form of the subcritical pitchfork bifurcation is

The fixed point structure in both supercritical and subcritical cases is shown in Fig. https://courses.physics.ucsd.edu/2011/Spring/physics221a/LECTURES/CH02_BIFURCATIONS.pdf

 

IV. Imperfect bifurcation

---

The imperfect bifurcation occurs when a symmetry-breaking term is added to the pitchfork. The normal form contains two control parameters:

Here, the constant h breaks the parity symmetry if \( u\,\mapsto -u . \)

This equation arises from a crude model of magnetization dynamics. Let M be the magnetization of a sample, and f(M) be the free energy. Assuming M is small, we can expand f(M) as

where H is the external magnetic field, and a and b are temperature-dependent constants. This is called the Landau expansion of the free energy. We assume b > 0 in order that the minimum of f(M) not lie at infinity. The dynamics of M(t) are modeled by with Γ > 0. Thus, the magnetization evolves toward a local minimum in the free energy. Note that the free energy is a decreasing function of time By rescaling \( M \equiv u\,M_0 \quad\mbox{with} \quad M_0 = (b\Gamma )^{-1/2} \) and defining \( k \equiv -a\Gamma \quad\mbox{and}\quad h \equiv \left( \Gamma^3 b \right)^{1/2} H , \) we obtain the normal form: Here, f(u) is a scaled version of the free energy.

Fixed points satisfy the equation

and correspond to extrema in f(u). By the fundamental theorem of algebra, this cubic polynomial may be uniquely factorized over the complex plane. Since the coefficients are real, the complex conjugate \( \overline{u} \) satisfies the same equation as u, hence there are two possibilities for the roots: either (i) all three roots are real, or (ii) one root is real and the other two are a complex conjugate pair. Clearly for k < 0 we are in situation (ii) since u³ - ku is then monotonically increasing for \( u \in \mathbb{R} , \) and therefore takes the value h precisely once for u real. For k > 0, there is a region \( h \in \left[ -h_c (k) , h_c (k) \right] \) over which there are three real roots. To find \( h_c (k) , \) we demand \( f'' (u) =0 \) as well as \( f' (u) =0 , \) which says that two roots have merged, forming an inflection point. One easily find \( h_c (k) = \frac{2}{3\sqrt{3}}\, k^{3/2} . \)

 

V. Landau Theory of Phase Transitions

---

 

VI. Single Neuron Equation

---

Consider a single neuron (nerve cell) that is receiving external input from surrounding cells in the brain, plus feedback from its own output,

Let x(t) denote the level of activity of the neuron at time t, normalized to be between 0 (low activity) and 1 (high activity). A simple model representing the change of activity level for such a cell is given by the differential equation

where E is the level of input from surrounding cells, θ is the cell's threshold, and S is called a response function. We will use a standard "digmoidal" response function This function increases from 0 and saturates at the value 1 as \( z \mapsto \infty . \)

The differential equation for x(t) becomes autonomous:

depending on parameter θ. The slope function is the difference of sigmoidal function S(z) and x. Since S(z) is always between 0 and 1, the slope function becomes negative for x > 1 and positive for x < 0. This means that any equilibrium solution \( x = S(x+E-\theta ) \) must lie between 0 and 1, and trajectories on the phase plane are always pointed down in half plane x > 1 and up in half plane x < 0.

The equilibrium solutions of this slope functions are constant values of x, where \( x = S(x+E-\theta ) . \) The figure above shows graphs of y = x and the response function \( y = S(x+E-\theta ) \) for values θ = 0.4, 0.7, and 1.0 with a = 10 and E=0.2. It is seen that for small θ there will be one equilibrium solution near x = 1 and for large θ there will be one stationary solution near x = 0. This seems reasonable because a high threshold means it takes a lot of input to produce very much activity. For θ in a middle range, however, there can be three equilibrium solutions. WE want to determine the two bifurcation values of θ at which the number of stationary solutions changes from 1 to 3 and from 3 to 1. First, we plot direction field: