Computer and Information Science » Numerical Analysis and Scientific Computing » "Numerical Simulation - From Theory to Industry", book edited by Mykhaylo Andriychuk, ISBN 978-953-51-0749-1, Published: September 19, 2012 under CC BY 3.0 license

# Bifurcation Analysis and Its Applications

By Canan Çelik Karaaslanlı
DOI: 10.5772/50075

Article top

## Overview

Figure 1. Bifurcation diagram corresponding to the saddle-node bifurcation

Figure 2. Bifurcation diagram corresponding to the transcritical bifurcation

Figure 3. Bifurcation diagram corresponding to the pitchfork bifurcation

Figure 4. Hopf bifurcation diagram

Figure 5. Bifurcation diagram corresponding to Supercritical Hopf bifurcation

Figure 6. Supercritical Hopf bifurcation

Figure 7. Diagram for Subcritical Hopf bifurcation

Figure 8. Subcritical Hopf bifurcation

Figure 9. The trajectory of prey density versus time with the initial condition N0=50. The graph of solutions of the model (29) when Ϣ=1.5<Ϣ0, where the equilibrium point E* is asymptotically stable.

Figure 10. The trajectory of predator density versus time with the initial condition P0=25. The graph of solutions of the model (29) when Ϣ=1.5<Ϣ0, where the equilibrium point E* is asymptotically stable.

Figure 11. The trajectories of S(t) versus time with the initial condition S0=25. The graph of solutions of the model (29) when Ϣ=1.5<Ϣ0, where the equilibrium point E* is asymptotically stable.

Figure 12. The trajectory of N,P,S in three dimension with the initial condition (50,25,25) when Ϣ=1.5<Ϣ0 for which the equilibrium point E* is asymptotically stable.

Figure 13. The trajectory of N(t) when Ϣ=2>Ϣ0.

Figure 14. The trajectory of P(t) when Ϣ=2>Ϣ0.

Figure 15. The trajectory of S(t) when Ϣ=2>Ϣ0.

Figure 16. When Ϣ=2>Ϣ0 a stable periodic orbit bifurcates from the equilibrium point E*.

# Bifurcation Analysis and Its Applications

Canan Çelik Karaaslanlı

## 1. Introduction

Continuous dynamical systems that involve differential equations mostly contain parameters. It can happen that a slight variation in a parameter can have significant impact on the solution. The main questions of interest in this chapter are: How to continue equilibria and periodic orbits of dynamical systems with respect to a parameter? How to compute stability boundaries of equilibria and limit cycles in the parameter space? How to predict qualitative changes in system’s behavior (bifurcations) occurring at these equilibrium points? This chapter will also cover the classification of bifurcations in terms of equilibria and periodic orbits. Especially it will present the specific bifurcation called ” Hopf bifurcation” which refers to the development of periodic orbits from stable equilibrium point, as a bifurcation parameter crosses a critical value. Since the theory of bifurcation from equilibria based on center manifold reduction and Poincaré-Normal forms, the direction of bifurcations for the mathematical models will also be explained using this theory. Finally, by introducing several software packages and numerical methods this chapter will also cover the techniques to determine and continue in some control parameters all local bifurcations of periodic orbits of dynamical systems and relevant normal form computations combined with the center manifold theorem, including periodic normal forms for periodic orbits.

In general, in a dynamical system, a parameter is allowed to vary, then the differential system may change. An equilibrium can become unstable and a periodic solution may appear or a new stable equilibrium may appear making the previous equilibrium unstable. The value of parameter at which these changes occur is known as ”bifurcation value” and the parameter that is varied is known as the ”bifurcation parameter”.

In this chapter, we also discuss several types of bifurcations, saddle node, transcritical, pitchfork and Hopf bifurcation. Among these types, we especially focus on Hopf bifurcation. The first three types of bifurcation occur in scalar and in systems of differential equations. The fourth type called Hopf bifurcation does not occur in scalar differential equations because this type of bifurcation involves a change to a periodic solution. Scalar autonomous differential equations can not have periodic solutions. Hopf bifurcation occurs in systems of differential equations consisting of two or more equations. This type is also referred to as a ”Poincare-Andronov-Hopf bifurcation”.

For a given system of differential equations first we shall consider the stability and the local Hopf bifurcation. By using the Hopf bifurcation theorem we prove the occurrence of the Hopf bifurcation. And then, based on the normal form method and the center manifold reduction introduced by Hassard et al.,[10], we derive the formulae determining the direction, stability and the period of the bifurcating periodic solution at the critical value of the bifurcation parameter. To verify the theoretical analysis, numerical simulations for bifurcation analysis are given in this chapter. For references see [1]-[22].

We also introduce the Hopf bifurcation for continuous dynamical systems and state the Hopf bifurcation theorem for these models. As it is well known, Hopf bifurcations occur when a conjugated complex pair of eigenvalues crosses the boundary of stability. In the time-continuous case, a limit cycle bifurcates. It has an angular frequency which is given by the imaginary part of the crossing pair. In the discrete case, the bifurcating orbit is generally quasi-periodic, except that the argument of the crossing pair times an integer gives just 2Ϟ. If we consider an ordinary differential equation (ODE) that depends on one or more parameters Ϗ

x'=f(x,Ϗ),

where, for simplicity, we assume Ϗ to be the only parameter. There is the possibility that under variation of Ϗ nothing interesting happens to Equation (1). There is only a quantitatively different behavior. Let us define Equation (1) to be structurally stable in the case there are no qualitative changes occurring. However, the ODE (Ordinary Differential Equation) might change qualitatively. At that point, bifurcations will have occurred.

Many of the basic principles for one dimensional systems apply also for two-dimensional systems. Let us define a two-dimensional system

x'=f(x,y,Ϗ),
y'=g(x,y,Ϗ),

where biologically we mostly interpret x as prey or resource and y as predator or consumer. Equilibria can be found by taking the equations equal to zero, i.e.,

f(x,y,Ϗ)=0,g(x,y,Ϗ)=0.

We have three possibilities for the stability of an equilibrium. Next to the stable and unstable equilibrium, there is the saddle equilibrium. A two-dimensional stable equilibrium is attracting in two directions, while a two-dimensional unstable equilibrium is repelling in two directions. A saddle point is attracting in one direction and repelling in the other direction. In the less formal literature saddles are often considered just unstable equilibria. A second remark is that also the dynamics of the system around the equilibria can differ.

The attracting or repelling can occur via straight orbits (a node) or via spiralling orbits (a spiral or focus). Note that, it is not possible to have a saddle focus in two dimensions. It is possible though in three of higher dimensional systems.

To prove the existence of Hopf bifurcation, we first obtain the Hopf bifurcation theorem hypothesis, i.e., the existence of purely imaginary eigenvalues of the corresponding characteristic equation with respect to the parameter Ϗ and also we prove the transversality condition

dϙdϏ(Ϗ=Ϗ0)0

at the critical value Ϗ0 where Hopf bifurcation occurs. Then based on the normal form approach and the center manifold theory introduced by Hassard et. al,[10], we derive the formula for determining the properties of Hopf bifurcation of the model.

Finally in this chapter, to support these theoretical results, we illustrate them by numerical simulations. In numerical analysis, generally MATLAB solver packages are used to analyze the dynamics of nonlinear models. In these solvers, differential equation systems are simulated by difference equations. In this chapter, we also give some examples from biology (such as well known predator-prey models with time delay) with numerical simulations and by graphing the solutions in two or three dimensions, we illustrate the occurrence of periodic solutions. For some examples see references by Çelik, [4], [5] and [6].

## 2. Basic concepts of bifurcation analysis

As it is stated above, in dynamical systems, a bifurcation occurs when a small smooth change made to the parameter values (the bifurcation parameters) of a system causes a sudden ” qualitative" or topological change in its behavior. Generally, at a bifurcation, the local stability properties of equilibria, periodic orbits or other invariant sets changes. It has two types;

Local bifurcations, which can be analyzed entirely through changes in the local stability properties of equilibria, periodic orbits or other invariant sets as parameters cross through critical thresholds; and

Global bifurcations, which often occur when larger invariant sets of the system ”collide” with each other, or with equilibria of the system. They cannot be detected purely by a stability analysis of the equilibria (fixed or equilibrium points, see the next section).

### 2.1. Equilibrium points

In dynamical systems, only the solutions of linear systems may be found explicitly. The problem is that in general real life problems may only be modeled by nonlinear systems. The main idea is to approximate a nonlinear system by a linear one (around the equilibrium point). Of course, we do hope that the behavior of the solutions of the linear system will be the same as the nonlinear one. But this is not always true. Before the linear stability analysis, we give some basic definitions below.

Definition (Equilibrium Point): Consider a nonlinear differential equation

x'(t)=f(x(t),u(t)),

where f is a function mapping RnxR3Rn. A point x̣ is called an equilibrium point if there is a specific ụϓRm such that

f(x(t),u(t))=0n.

Suppose x̣ is an equilibrium point (with the input ụ). Consider the initial condition x(0)=x̣, and applying the input u(t)=ụ for all tt0, then resulting solution x(t) satisfies

x(t)=x̣,

for all tt0. That is why it is called an equilibrium point or solution.

Example: As an example, consider the logistic growth equation (the rate of population density)

x'=rx(1-xK),

where x(t) denotes the population density at time t, r and K are positive constants, K is the carrying capacity. Then by setting right hand side function equal to zero,

f(x)=rx(1-xK)=0,

we obtain two equilibrium points x=0 and x=K.

### 2.2. Linear stability analysis

Linear stability of dynamical equations can be analyzed in two parts: one for scalar equations and the other for two dimensional systems;

#### 2.2.1. Linear stability analysis for scalar equations

To analyze the ODE

x'=f(x)

locally about the equilibrium point x=x̣, we expand the function f(x) in a Taylor series about the equilibrium point x̣. To emphasize that we are doing a local analysis, it is customary to make a change of variables from the dependent variable x to a local variable. Now let

x(t)=x̣+ϓ(t),

where it is assumed that ϓ(t)1, so that we can justify dropping all terms of order two and higher in the expansion. Substituting x(t)=x̣+ϓ(t) into the RHS of the ODE yields;

f(x(t))=f(x̣+ϓ(t))=f(x̣)+f'(x̣)ϓ(t)+f''(x̣)ϓ2(t)2+...=0+f'(x̣)ϓ(t)+O(ϓ2),

and dropping higher order terms, we obtain

f(x)f'(x̣)ϓ(t).

Note that dropping these higher order terms is valid since ϓ(t)1. Now substituting x(t)=x̣+ϓ(t) into the LHS of the ODE,

ϓ'(t)=f'(x̣)ϓ(t).

The goal is to determine if we have growing or decaying solutions. If the solutions grows, then the equilibrium point is unstable. If the solution decays, then the fixed point is stable.

To determine whether or not the solution is stable or unstable we simply solve the ODE and get the solution as

ϓ(t)=ϓ0exp(f'(x̣)t),

where ϓ0 is a constant. Hence, the solution is growing if f'(x̣)>0 and decaying if f'(x̣)<0. As a result, the equilibrium point is stable if f'(x̣)<0, unstable if f'(x̣)>0 as it is stated in the following theorem.

Theorem: Suppose for scalar differential equation

x'=f(x),

the derivative function f' is continuous on an open interval I where the equilibrium point x̣ϓI. Then the equilibrium point x̣ is locally stable if f'(x̣)<0 and it is unstable if f'(x̣)>0.

If the equilibrium point is stable and in addition

limtx(t)=x̣,

then it is called asymptotically stable equilibrium point.

Example: Determine the stability of the fixed points to the Logistic growth equation

N'=fN=rN1-NK,wherer>0.

We first find the equilibrium points by setting

fṆ=0,whichyieldstwopointsṆ=0,K.

Next we compute f'(N) and evaluate it at the equilibrium points.

f'(N)=r(1-2NK)atṆ=0,f'(0)=r>0,soitisunstableatṆ=K,f'(K)=-r<0,soitislocallystable.

#### 2.2.2. Linear stability analysis for systems

Consider the two dimensional nonlinear system

x'=f(x,y),y'=g(x,y),

and suppose that (x̣,ỵ) is a steady state (equilibrium point), i.e.,

fx̣,ỵ=0andg(x̣,ỵ)=0.

Now let’s consider a small perturbation from the steady state (x̣,ỵ)

x=x̣+u,y=ỵ+v,

where u and v are understood to be small as u1 and v1. Is is natural to ask whether u and v are growing or decaying so that x and y will move away form the steady state or move towards the steady states. If it moves away, it is called unstable equilibrium point, if it moves towards the equilibrium point, then it is called stable equilibrium point.As in scalar equations, by expanding the Taylor’s series for f(x,y) and g(x,y);

u'=x'=f(x,y)=f(x̣+u,ỵ+v)=fx̣,ỵ+fxx̣,ỵu+fyx̣,ỵv+higherorderterms...=fxx̣,ỵu+fyx̣,ỵv+higherorderterms....

Similarly,

v'=y'=g(x,y)=g(x̣+u,ỵ+v)=gx̣,ỵ+gxx̣,ỵu+gyx̣,ỵv+higherorderterms...=gxx̣,ỵu+gyx̣,ỵv+higherorderterms....

Since u and v are assumed to be small, the higher order terms are extremely small, we can neglect the higher order terms and obtain the following linear system of equations governingthe evolution of the perturbations u and v,

u'v'=fx(x̣,ỵ)fy(x̣,ỵ)gx(x̣,ỵ)gy(x̣,ỵ)uv,

where the matrix fxfygxgy is called Jacobian matrix J of the nonlinear system.The above linear system for u and v has the trivial steady state (u,v)=(0,0), and the stability of this trivial steady state is determined by the eigenvalues of the Jacobian matrix, as follows:

Theorem: An equilibrium point (x̣,ỵ) of the differential equation is stable if all the eigenvalues of J, the Jacobian evaluated at (x̣,ỵ) have negative real parts. The equilibrium point is unstable if at least one of the eigenvalues has a positive real part.

As a summary,

Asymptotically stable : A critical point is asymptotically stable if all eigenvalues of the jacobian matrix J are negative, or have negative real parts.

Unstable: A critical point is unstable if at least one eigenvalue of the jacobian matrix J is positive, or has positive real part.

Stable (or neutrally stable) : Each trajectory move about the critical point within a finite range of distance.

Definition(Hyperbolic point): The equilibrium is said to be hyperbolic if all eigenvalues of the jacobian matrix have non-zero real parts.

Hyperbolic equilibria are robust(i.e., the system is structurally stable): Small perturbations of order do not change qualitatively the phase portrait near the equilibria. Moreover, local phase portrait of a hyperbolic equilibrium of a nonlinear system is equivalent to that of its linearization. This statement has a mathematically precise form known as the Hartman-Grobman. This theorem guarantees that the stability of the steady state (x̣,ỵ) of the nonlinear system is the same as the stability of the trivial steady state (0,0) of the linearized system.

Definition(Non-Hyperbolic point): If at least one eigenvalue of the Jacobian matrix is zero or has a zero real part, then the equilibrium is said to be non-hyperbolic.

Non-hyperbolic equilibria are not robust (i.e., the system is not structurally stable). Small perturbations can result in a local bifurcation of a non-hyperbolic equilibrium, i.e., it can change stability, disappear, or split into many equilibria. Some refer to such an equilibrium by the name of the bifurcation (See the section below).

Example: Consider the following nonlinear autonomous system

x(t)=y(t)[x(t)-y(t)],
y(t)=x(t)[2-y(t)].

The equilibria are the points (x̣,ỵ)=(0,0) and (x̣,ỵ)=(1,2) and the Jacobian matrix is

J=yx-12-y-x.

We compute the Jacobian at the equilibrium point (0,0) where J(0,0)=0-120 which implies that the eigenvalues are purely imaginary

ϙ1.2=±2i,

by solving the characteristic equation

det(J-ϙI)=0.

Since the point is non-hyperbolic, the linearized system can not tell about the stability. Later on we will show that this is a center.

For the equilibrium point (1,2), the Jacobian J(1,2)=200-1 thus the point is locally unstable (as

ϙ1=1 and ϙ2=-1,

where one of the eigenvalues is strictly positive). Since it is a hyperbolic equilibrium point, the stability of fixed point is the same as in the linearized system. So it is also unstable.

## 3. Bifurcation analysis

Consider a family of ODE’s that depend on one parameter ϙ

x'=f(x,ϙ),

where f:Rn+1Rn is analytic for ϙϓR, xϓRn. Let x=x0(ϙ) be a family of equilibrium points of equation (5), i.e., f(x0(ϙ),ϙ)=0. Now let’s set

z=x-x0(ϙ).

Then

z'=A(ϙ)z+O(|z|2),

where A(ϙ)=fx(x0(ϙ),ϙ).

Let ϙ1,ϙ2,...ϙn(ϙ) be the eigenvalues of A(ϙ). If, for some i, Reϙi(ϙ) changes sign at ϙ=ϙ0, we say that ϙ0 is a bifurcation point of equation (5).

### 3.1. Bifurcation in one dimension

We may assume that n=1 so that f:R2R1, and x0(ϙ) is a real valued analytic function of ϙ provided

ϙ1(ϙ)=fx(x0(ϙ),ϙ)=A(ϙ)=0.

Therefore, the equilibrium point is asymptotically stable if ϙ1(ϙ)<0, and unstable if ϙ1(ϙ)>0. This implies that ϙ0 is a bifurcation point if ϙ1(ϙ0)=0. Hence, bifurcation points (x0(ϙ),ϙ)are solutions of

f(x,ϙ)=0 andfx(x,ϙ)=0.

The most common bifurcation types are illustrated by the following examples.

A saddle-node bifurcation or tangent bifurcation is a collision and disappearance of two equilibria in dynamical systems. In autonomous systems, this occurs when the critical equilibrium has one zero eigenvalue. This phenomenon is also called fold or limit point bifurcation.

Now consider the dynamical system defined by

x'=a-x2, for a isreal.

An equilibrium solution (where x'=0) is simply x=±a. Therefore,

if a<0, then we have no real solutions,if a>0, then we have two real solutions.

We now consider each of the two solutions for a>0, and examine their linear stability in the usual way. First, we add a small perturbation:

x=x̣+ϓ.

Substituting this into the equation yields

dϓdt=(a-x̣2)-2x̣ϓ-ϓ2,

and since the term in brackets on the RHS is trivially zero, therefore

dϓdt=-2x̣ϓ,

which has the solution

ϓ(t)=Aexp(-2x̣t).

From this, we see that for x=+a, |x|0 as t (linear stability); for x=-a,|x|0 as t (linear stability).

As sketched in the “ bifurcation diagram” below, therefore, the saddle node bifurcation at a=0 corresponds to the creation of two new solution branches. One of these is linearly stable, the other is linearly unstable.

### Figure 1.

Bifurcation diagram corresponding to the saddle-node bifurcation

#### 3.1.2. Transcritical bifurcation

In a transcritical bifurcation, two families of fixed points collide and exchange their stability properties. The family that was stable before the bifurcation is unstable after it. The other fixed point goes from being unstable to being stable.

Now consider the dynamical system

dxdt=az-bx2, for x,a,b real.

Again, a and b are control parameters. We can find two steady states (x'=0) to this system

x=x̣1=0,a,bx=x̣2=ab,a,b,b0.

We now examine the linear stability of each of these states in turn, following the usual procedure.

For the state x̣1, we add a small perturbation

x=x̣1+ϓ,

which yields

dϓdt=aϓ-bϓ2

with the linearized form

dϓdt=aϓ

has the solution

e(t)=Aexp(at).

Therefore, perturbations grow for a > 0 and decay for a < 0. So

thestatex̣1isunstableifa>0,thestatex̣1isstableifa<0.

Now for the state x̣2, we add a small perturbation

x=x̣2+ϓ,

which yields the linearized form

dϓdt=-aϓ

has the solution

e(t)=Aexp(-at).

Therefore, perturbations grow for a > 0 and decay for a < 0. So

thestatex̣2=abislinearlystableifa>0,thestatex̣2=abislinearlyunstableifa<0.

It can be easily seen that the bifurcation point a=0 corresponds to an exchange of stabilities between the two solution branches.

### Figure 2.

Bifurcation diagram corresponding to the transcritical bifurcation

#### 3.1.3. The pitchfork bifurcation

In pitchfork bifurcation one family of fixed points transfers its stability properties to two families after or before the bifurcation point. If this occurs after the bifurcation point, then pitchfork bifurcation is called supercritical. Similarly, a pitchfork bifurcation is called subcritical if the nontrivial fixed points occur for values of the parameter lower than the bifurcation value. In other words, the cases in which the emerging nontrivial equilibria are stable are called supercritical whereas the cases in which these equilibria are called subcritical.

Consider the dynamical system

x'=ax-bx3,for a,b real.

As usual, a and b are external control parameters. Steady states, for which x'=0 are as follows:

x=x̣1=0,a,b,x=x̣2=-a/b,fora/b>0,x=x̣3= a/b,fora/b>0.

Note that the equilibrium points x̣2 and x̣3 only exist when a>0 if b>0 and for a<0 if b<0.

As usual, we now examine the linear stability of each of these steady states in turn. (This can be done for a general b). First we write the perturbation for x̣1=0,

x=x̣1+ϓ

that yields the linearized equation

dϓdt=aϓ,

with the solution

e(t)=Aexp(at).

So we see that

thestatex̣1=0islinearlyunstableifa>0,thestatex̣2=0islinearlystableifa<0.

For the states x=x̣2 and x=x̣3, setting

x̣=±a/b+ϓ,dϓdt=aϓ-3bx̣2ϓ,

with the solution

e(t)=Aexp(ct)wherec=-2a.

Thus it is obvious that

thestatesx̣2andx̣3arelinearlystableifa>0,thestatesx̣2andx̣3arelinearlyunstableifa<0,

#### Figure 3.

Bifurcation diagram corresponding to the pitchfork bifurcation

### 3.2. Bifurcation in two dimension

#### 3.2.1. Hopf bifurcation

Definition: A Hopf or Poincare-Andronov-Hopf bifurcation is a local bifurcation in which a fixed point of a dynamical system loses stability as a pair of complex conjugate eigenvalues of linearization around the fixed point cross the imaginary axis of the complex plane.

#### Figure 4.

Hopf bifurcation diagram

### 3.3. Hopf bifurcation theorem

Consider the two dimensional system

dxdt=fx,y,Ϣ,
dydt=gx,y,Ϣ,

where Ϣ is the parameter and suppose that x¯Ϣ,y¯Ϣ is the equilibrium point and ϏϢ±iϐϢ are the eigenvalues of the Jacobian matrix which is evaluated at the equilibrium point.

In addition let’s assume that the change in the stability of the equilibrium point occurs at Ϣ=Ϣ* where ϏϢ*=0.

First the system is transformed so that the equilibrium is at the origin and the parameter Ϣ at Ϣ*=0 gives purely imaginary eigenvalues. System (7) is rewritten as follows;

dxdt=a11Ϣx+a12Ϣy+f1(x,y,Ϣ),
dydt=a21Ϣx+a22Ϣy+g1x,y,Ϣ.

The linearization of the system (7) about the origin is given by dXdt=J(Ϣ)X, where X=xy and

J(Ϣ)=a11Ϣa12Ϣa21Ϣa22Ϣ

is the Jacobian matrix evaluated at origin.

Theorem (Hopf bifurcation theorem)

Let f1 and g1, in system (8) have continuous third order partial derivatives in x and y. Suppose that the origin is an equilibrium point of (8) and that the Jacobian matrix JϢ as above, is valid for all sufficiently small |Ϣ|. Moreover, assume that the eigenvalues of matrix JϢ are ϏϢ±iϐϢ where Ϗ0=0,ϐ00 such that the eigenvalues cross the imaginary axis with nonzero speed, i.e., dϏdϢ|Ϣ=00.

Then in any open set U containing the origin in R2 and for any Ϣ0>0, there exists a value Ϣ̣,|Ϣ̣|<Ϣ0 such that the system of differential equations (8) has a periodic solution for Ϣ=Ϣ̣ in U. (Allen, L.J.S).

Note: The Hopf bifurcation requires at least a two dimensional differential equation system to appear.

Definition: The bifurcation stated in the Hopf bifurcation theorem is called ” supercritical” if the equilibrium point (0,0) is asymptotically stable when Ϣ=0 (at the bifurcation point) and it is called ” subcritical” if the equilibrium point (0,0) is negatively asymptotically stable (as t-) when Ϣ=0.

In a supercritical Hopf bifurcation, the limit cycle grows out of the equilibrium point. In other words, right at the parameters of the Hopf bifurcation, the limit cycle has zero amplitude, and this amplitude grows as the parameters move further into the limit-cycle. (See the figure below)

### Figure 5.

Bifurcation diagram corresponding to Supercritical Hopf bifurcation

### Figure 6.

Supercritical Hopf bifurcation

However in a subcritical Hopf bifurcation, there is an unstable limit cycle surrounding the equilibrium point, and a stable limit cycle surrounding that. The unstable limit cycle shrinks down to the equilibrium point, which becomes unstable in the process. For systems started near the equilibrium point, the result is a sudden change in behavior from approach to a stable focus, to large-amplitude oscillations.(See the figure below).

### Figure 7.

Diagram for Subcritical Hopf bifurcation

### Figure 8.

Subcritical Hopf bifurcation

Example: Consider the two dimensional system

dxdt=x2-Ϗ,
dydt=-(x2+1)y,

where Ϗ is a parameter. When we compute the equilibrium points depending on parameter Ϗ;

If Ϗ<0, then there is no x-nullclines, hence the system has no equilibrium points.If Ϗ=0, then the system has exactly one equilibrium point at (0,0).If Ϗ>0, then the system has two equilibrium points (-Ϗ,0) and (Ϗ,0).

Then the Jacobian matrix is

J=2x0-2xy-(x2+1),

where at the equilibrium points J(0.0)=000-1,J(Ϗ,0)=-2Ϗ00-Ϗ-1 and J(Ϗ,0)=2Ϗ00-Ϗ-1.

For Ϗ=0, there will be a line equilibrium (since one of the eigenvalues is zero) and for Ϗ>0, the point (-Ϗ,0) is a sink and (Ϗ,0) is a saddle point so that Ϗ=0 is the bifurcation point for this differential equation system.

Example: Consider the two dimensional system

dxdt=y+Ϗx,
dydt=-x+Ϗy,

where Ϗ is the bifurcation parameter. We can easily show that the conditions of the Hopf Bifurcation theorem hold. In this system f1 and f2 are zero. Then the Jacobian matrix is

J=Ϗ1-1Ϗ

for which the eigenvalues are ϙ1,2=Ϗ±i where Reϙ(Ϗ)=Ϗ and the imaginary part ImϙϏ=±1.It follows that Reϙ(0)=0 and Imϙ(0)0. and also

dReϙ(Ϗ)dϏ|Ϗ=0=10.

Hence, we conclude that there exists a periodic solution for Ϗ=0 in every neighborhood of origin.

## 4. Center Manifold teorem

Let f(0)=0, for the dynamical system
x'=f(x),xRn,

and let the eigenvalues of the Jacobian matrix be ϙ1,ϙ2,...,ϙn. Suppose that, the real parts of the eigenvalues are zero and if not, suppose there are n+ numbers of eigenvalues with Reϙ>0,n0 number of eigenvalues with Reϙ=0 and n- number of eigenvalues with Reϙ<0. Let Tc be the eigenspace on imaginary axis corresponding to n0 eigenvalues. The eigenvalues on the imaginary axis Reϙ=0 are called the critical eigenvalues as on the eigenspace Tc. And suppose the function Ϥt denote the flow corresponding to the equation (11).

With these assumption, we state the Center Manifold theorem as follows;

Theorem: (Center Manifold theorem)

There exists a locally invariant C center manifold Wlocc0 such that

Wlocc0={(x,y):y=h(x);|x|<ϒ,h(0)=0;DJ(0)=0}

such that the dynamics of the system

x'=Acx+r1(x,y),y'=Asy+r2(x,y),

(where Acand As are are the blocks in the canonical form whose diagonals contain the eigenvalues with Reϙ=0 and Reϙ<0; respectively) restricted to the center manifold are given by

x'=Ac(x)+r1(x,h(x)).

And the manifold Wlocc is called center manifold.

Remark: Center manifolds are not unique.

### 4.1. Center Manifold reduction for two dimensional systems

Example: Consider the two dimensional system of differential equations

x'=xy,
y'=-y+x2.

The only equilibrium point is (0,0), we linearize around that and obtain

J(0,0)=000-1.

Now we look for y=h(x)=ax2+bx3+cx4+dx5+O(x6). Then,

y'=h'(x)x'=xh'(x)h(x)=2a2x4+5abx5+O(x6)

and on the other hand,

y'=-h(x)-x2=-(a+1)x2-bx3-cx4-dx5+O(x6).

Comparing the two expressions, we deduce that a=-1,b=0,c=-2,d=0 and the center manifold reduction takes the form

y=h(x)=-x2-2x4+O(x7).

Hence for the last equation x=0 is asymptotically stable and therefore (0,0) is asymptotically stable for the original system.

Example: Consider the two dimensional system

x'=x2y-x5,
y'=-y+x2.

Again (0,0) is an equilibrium point and the jacobian matrix for the linearized system is

J(0,0)=000-1.

Consider the transformation y=h(x)=ax2+bx3+cx4+dx5+O(x6), which leads

y'=h'(x)x'=2a2x5+[2a(b-1)+3ab]x6+O(x7),and y'=-h(x)+x2=-(a+1)x2-bx3+O(x4)=2a2x4+5a-bxbx5+O(x6),

from which we deduce that a=1 and b=0 and

x'=x4+O(x5).

For this reduced equation, x=0 is unstable and hence, (0,0) is also unstable for the original system.

### 4.2. Center Manifold reduction for Hopf bifurcation

The aim of this section is to give a formal framework for the analytical bifurcation analysis of Hopf bifurcations in delay differential equations

x'=f(x,Ϣ),xϓR3

with a single fixed time delay Ϣ to be chosen as a bifurcation parameter. Characteristic equations of the delay differential equation form (14) are often studied in order to understand changes in the local stability of equilibria of certain delay differential equations. It is therefore important to determine the values of the delay at which there are roots with zero part. We give a general formalization of these calculations and determine closed form algebraic equations where the stability and amplitude of periodic solutions close to bifurcation can be calculated.

We shall determine the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions by applying the normal form theory and the center manifold theorem by Hassard et al.,[10], and throughout this section, we assume that the three dimensional system of delay differential equations (14) undergoes Hopf bifurcations at the positive equilibrium (N0*,P0*,S0*) at Ϣ=Ϣk, and iϧ1 is the corresponding purely imaginary root of the characteristic equation at the positive equilibrium (N0*,P0*,S0*). For the sake of simplicity, we use the notation iϧ for iϧ1.

We first consider the system (14) by the transformation

x1=N-N0*,x2=P-P0*,x3=S-S0*,t=tϢ,Ϣ=Ϣk+Ϛ

which is equivalent to the following Functional Differential Equation(FDE) system in C=C([-1,0],R3),

x̥(t)=LϚ(xt)+f(Ϛ,xt),

where x(t)=(x1(t),x2(t),x3(t))TϓR3, and LϚ:CR3,f:R×CR3 are given respectively, by

LϚ(ϳ)=(Ϣk+Ϛ)a1a2a3a4a5a6a7a8a9ϳ1(0)ϳ2(0)ϳ3(0)+(Ϣk+Ϛ)b1b2b3b4b5b6b7b8b9ϳ1(-1)ϳ2(-1)ϳ3(-1),andf(Ϛ,ϳ)=(Ϣk+Ϛ)f11f12f13.

By Riesz representation theorem, there exists a function ϕ(ϖ,Ϛ) of bounded variation for ϖϓ-1,0], such that

LϚϳ=-10dϕ(ϖ,0)ϳ(ϖ),forϳϓC.

Indeed we may take

ϕ(ϖ,Ϛ)=(Ϣk+Ϛ)a1a2a3a4a5a6a7a8a9ϒ(ϖ)-(Ϣk+Ϛ)b1b2b3b4b5b6b7b8b9ϒ(ϖ+1),

where ϒ is the Dirac delta function. For ϳϓC1([-1,0],R3), define

A(Ϛ)ϳ=(ϖ)dϖϖϓ-1,0),-10dϕ(Ϛ,s)ϳ(s) ϖ=0andR(Ϛ)ϳ=0 ϖϓ-1,0),f(Ϛ,ϳ),ϖ=0.

Then the system (15) is equivalent to

x̥t=A(Ϛ)xt+R(Ϛ)xt,

where xt(ϖ)=x(t+ϖ) for ϖϓ-1,0). For ϦϓC1([-1,0],(R3)*), define

A*Ϧ(s)=-dϦ(s)dssϓ(0,1],-10dϕT(t,0)Ϧ(-t) s=0

and a bilinear inner product

Ϧ(s),ϳ(ϖ)=Ϧ̣(0)ϳ(0)--10Ϝ=0ϖϦ̣(Ϝ-ϖ)dϕ(ϖ)ϳ(Ϝ)dϜ,

where ϕ(ϖ)=ϕ(ϖ,0). Then A(0) and A* are adjoint operators. Suppose that q(ϖ) and q*(s) are eigenvectors of A and A* corresponding to iϧϢk and -iϧϢk, respectively. Then suppose that q(ϖ)=(1,ϐ,ϑ)TeiϧϢkϖ is the eigenvector of A(0) corresponding to iϧϢk, then A(0)q(ϖ)=iϧϢkq(ϖ). Then in the following, we use the theory by Hassard et al.,[10], to compute the coordinates describing center manifold C0 at Ϛ=0. Define

z(t)=q*,xt, W(t,ϖ)=xt-2Rez(t)q(ϖ).

On the center manifold C0, we have

W(t,ϖ)=W(z(t),ẓ(t),ϖ)=W20(ϖ)z22+W11(ϖ)zẓ+W02(ϖ)ẓ22+...,

where z and ẓ are local coordinates for center manifold C0 in the direction of q and q̣*. Note that W is real if xt is real. We consider only real solutions. For the solution xtϓC0, since Ϛ=0 and (15), we have

z̥=iϧϢkz+q*(ϖ),F(0,W(z,ẓ,ϖ)+2Rezq(ϖ))=iϧϢkz+q̣*(0)F(0,W(z,ẓ,0)+2Rezq(0))=defiϧϢkz+q̣*(0)F0(z,ẓ)=iϧϢkz+g(z,ẓ),

where

g(z,ẓ)=q̣*(0)F0(z,ẓ)=g20z22+g11zẓ+g02ẓ22+g21z2ẓ2+....

By using (17), we have xt(x1t(ϖ),x2t(ϖ),x3t(ϖ))=W(t,ϖ)+zq(ϖ)+zq¯(ϖ), q(ϖ)=(1,ϐ,ϑ)TeiϧϢkϖ, and

x1t(0)=z+ẓ+W20(1)(0)z22+W11(1)(0)zẓ+W02(1)(0)ẓ22+O((z,ẓ)3),x2t(0)=ϐz+ϐ̣z¯+W20(2)(0)z22+W11(2)(0)zẓ+W02(2)(0)ẓ22+O((z,ẓ)3),x3t(0)=ϑz+ϑ¯Ϗ¯+W20(3)(0)z22+W11(3)(0)zẓ+W02(3)(0)ẓ22+O((z,ẓ)3),x1t(-1)=ze-iϧϢkϖ+ẓeiϧϢkϖ+W20(1)(-1)z22+W11(1)(-1)zẓ+W02(1)(-1)ẓ22+O((z,ẓ)3).

From the definition of F(Ϛ,xt), we have

g(z,ẓ)=q̣*(0)f0(z,ẓ)=ḌϢk(1,ϐ̣*,ϑ̣*)f110f120f130

and we evaluate g(z,ẓ).

To determine g21, we need to compute W20(ϖ) and W11(ϖ). By (15) and (18), we have

W̥=x̥t-z̥q+zq¯.
={AW-2Re{q̣*(0)f0q(ϖ)}, ϖϓ-1,0)AW-2Re{q̣*(0)f0q(ϖ)}+f0ϖ=0, =defAW+H(z.ẓ,ϖ),

where

H(z.ẓ,ϖ)=H20(ϖ)z22+H11(ϖ)zẓ+H02(ϖ)ẓ22+....

Note that on the center manifold C0 near to the origin,

W̥=Wzz̥+Wẓẓ.

Thus we obtain,

(A-2iϧϢk)W20(ϖ)=-H20(ϖ), AW11(ϖ)=-H11(ϖ).

By using (19), for ϖϓ-1,0),

H(z.ẓ,ϖ)=q̣*(0)f0q(ϖ)-q*(0)f0(0)q̣(ϖ)=-gq(ϖ)-g̣q̣(ϖ).

Comparing the coefficients with (20), we obtain the following

H20(ϖ)=-g20q(ϖ)-g̣02q̣(ϖ), H11(ϖ)=-g11q(ϖ)-g̣11q̣(ϖ).

From (22) and (24) and the definition of A, we get

W̥20(ϖ)=2iϧϢkW20(ϖ)-g20q(ϖ)-g̣02q̣(ϖ).

Noticing q(ϖ)=q(0)eiϧϢkϖ, we evaluate W20(ϖ) by E1=(E1(1),E1(2),E1(3))ϓR3 is a constant vector. From the definition of A and (22), we obtain

-10dϕ(ϖ)W20(ϖ)=2iϧϢkW20(0)-H20(0),
and
-10dϕ(ϖ)W11(ϖ)=-H11(0),

where dϕ(ϖ)=ϕ(ϖ,0). Next we compute W20(ϖ) and W11(ϖ) from (25) and (26) and determine the following values to investigate the qualities of bifurcating periodic solution in the center manifold at the critical value Ϣk. For this purpose, we express the direction of Hopf bifurcation in terms of gij's and eigenvalues ϙϢk. And then we can evaluate the following values;

c1(0)=i2ϧϢk(g20g11-2|g11|2-|g02|23)+g212,Ϛ2=-Re{c1(0)}Re{ϙ'(Ϣk)},
ϐ2=2Re{c1(0)},T2=-Im{c1(0)}+Ϛ2Im{ϙ'(Ϣk)}ϧϢk

and we state this as in the following theorem.

Theorem: Ϛ2 determines the direction of Hopf bifurcation; if Ϛ2>0, then the Hopf bifurcation is supercritical and the bifurcating periodic solutions exist for Ϣ>Ϣ0, if Ϛ2<0, then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exist for Ϣ<Ϣ0. ϐ2 determines the stability of the bifurcating periodic solutions; bifurcating periodic solutions are stable if ϐ2<0, unstable if ϐ2>0. T2 determines the period of the bifurcating solution; the period increases if T2>0, it decreases if T2<0.

In the following section, we shall give a numerical example to verify the theoretical results.

#### 4.2.1. Numerical example of Center Manifold reduction

Consider the following system with discrete time delay Ϣ;

dN(t)dt=r1N(t)-ϓP(t)N(t),dP(t)dt=P(t)(r2-ϖS(t)N(t-Ϣ)),dS(t)dt=ϏP(t)-ϏS(t).

Now we present some numerical simulations by using MATLAB(7.6.0) programming (Çelik-3). We simulate the predator-prey system (28) by choosing the parameters r1=0.45, r2=0.1, ϖ=0.05, ϓ=0.03, and Ϗ=1, i.e., we consider the following system

dN(t)dt=0.45N(t)-0.03P(t)N(t),dP(t)dt=P(t)(0.1-0.05S(t)N(t-Ϣ)),dS(t)dt=P(t)-S(t),

which has only one positive equilibrium E*=(N0*,P0*,S0*)=(7.5,15,15). By algorithms in the previous sections, we obtain Ϣ0=1.5663, ϧ1=0.0045,z1=ϧ12 and g'(z1)=1.1944×10-4>0 which leads to dReϙ(Ϣ0)dϢ=5.8982>0. So by the theorem above, the equilibrium point E* is asymptotically stable when Ϣϓ[0,Ϣ0)=[0,1.5663) and unstable when Ϣ>1.5663 and also Hopf bifurcation occurs at Ϣ=Ϣ0=1.5663 as it is illustrated by computer simulations.

By the theory of Hassard et al.,[10], as it is discussed in previous section, we also determine the direction of Hopf bifurcation and the other properties of bifurcating periodic solutions. From the formulae in Section 5.2 we evaluate the values of Ϛ2, ϐ2 and T2 as

Ϛ2=0.0981>0, ϐ2=-0.5785<0, T2=7.1165>0,

from which we conclude that Hopf bifurcation of system (29) occurring at Ϣ0=1.5663 is supercritical and the bifurcating periodic solution exists when Ϣ crosses Ϣ0to the right, and also the bifurcating periodic solution is stable.

In computer simulations, the initial conditions are taken as (N0,P0,S0)=(50,25,25) and MATLAB DDE (Delay Differential Equations) solver is used to simulate the system (29). We first take Ϣ=1.5<Ϣ0 and plot the density functions N(t), P(t) and S(t) in Figs.9,10,11 respectively which shows the positive equilibrium is asymptotically stable for Ϣ<Ϣ0.. Moreover in Fig.12, we illustrate the asymptotic stability in three dimension.

### Figure 9.

The trajectory of prey density versus time with the initial condition N0=50. The graph of solutions of the model (29) when Ϣ=1.5<Ϣ0, where the equilibrium point E* is asymptotically stable.

### Figure 10.

The trajectory of predator density versus time with the initial condition P0=25. The graph of solutions of the model (29) when Ϣ=1.5<Ϣ0, where the equilibrium point E* is asymptotically stable.

### Figure 11.

The trajectories of S(t) versus time with the initial condition S0=25. The graph of solutions of the model (29) when Ϣ=1.5<Ϣ0, where the equilibrium point E* is asymptotically stable.

### Figure 12.

The trajectory of N,P,S in three dimension with the initial condition (50,25,25) when Ϣ=1.5<Ϣ0 for which the equilibrium point E* is asymptotically stable.

### Figure 13.

The trajectory of N(t) when Ϣ=2>Ϣ0.

### Figure 14.

The trajectory of P(t) when Ϣ=2>Ϣ0.

### Figure 15.

The trajectory of S(t) when Ϣ=2>Ϣ0.

### Figure 16.

When Ϣ=2>Ϣ0 a stable periodic orbit bifurcates from the equilibrium point E*.

However in Figs.13,14,15 and 16 below, we take Ϣ=2>Ϣ0 sufficiently close to Ϣ0 which illustrates the existence of bifurcating periodic solutions from the equilibrium point E*.

### 4.3. Exercises

1. Consider the Van der Pol equation

x̦-(1-x2)x̥+x=0

and convert this into a system and check the stability of fixed points.

2. Let y0 be an equilibrium point of the equation

dydt=f(y).Let g(y) be a positive function. Deduce the stability of y0 as an equilibrium solution of the equationdydt=f(y)g(y)

from its stability as an equilibrium solution of the equation dydt=f(y). Repeat the same question if g(y)is a negative function.

3. For the following nonlinear system

x'=-x3-y,y'=x-y3,

determine the stability of fixed points.

4. For the following nonlinear system

x'=y-xy2,y'=-x+x2y,

classify the fixed points.

5. Analyze the bifurcation properties of the following problems choosing r as bifurcation parameter,

1. x'=-x+ϐtanhx,
2. x'=rx-4x3,
3. x'=rx-sin(x),
4. x'=rx+4x3,
5. x'=rx-sinh(x),
6. x'=x+rx1+x2.

6. Find the equilibrium points and identify the bifurcation in the following system, and sketch the appropriate bifurcation diagram and phase portraits:

dxdt=(1+Ϣ)y-x-2x3,dydt=x-y-y3.

Then compute the extended center manifold near the bifurcation point by choosing Ϣ as bifurcation parameter.

7. Show that the following system is structurally unstable,

dxdt=x+y-x2+y2,dydt=-2x-y+xy,

and the following system is structural stable. Explain our reason.

dxdt=2x+y-x(x2+y2),dydt=-x+2y-y(x2+y2).

8.

dxdt=xy,dydt=1-y2-x2,

for the above system, classify the fixed points.

9. For rR, consider the differential equation

x'=rx-2x2+x3

on the real line.

1. Show that x*=0 is a fixed point for any value of the parameter r, and determine its stability. Hence identify a bifurcation point r1.

2. Show that for certain values of the parameter r there are additional fixed points. For which values of r do these fixed points exist? Determine their stability and identify a further bifurcation points r2.

3. Using a Taylor expansion of the differential equation above, determine the normal form of the bifurcation at r1. What type of bifurcation takes place.

4. Similarly, determine the normal from of the bifurcation at r2. What type bifurcation takes place?

5. Sketch the bifurcation diagram for all values of r and x*. (Use a full line to denote a curve of stable fixed points, and a dashed line for a curve of unstable fixed points).

10. Consider the system of differential equations

x'=y-xy2,y'=-x+yx2,

in the (x,y) plane.

1. Determine all fixed points of the system.

2. Let r2=x2+y2. Show that r' =0 Considering the result of (c), what does this imply for the trajectories?
3. Sketch the phase portrait in the (x,y) plane, including trajectories through (1,0) and (2,0). Whichfixed point does the trajectory through (2,0) approach?

11. Consider the Lorenz system (the model of heat convection by Rayleigh-Benard occurring in the earth’s atmosphere) in three dimension as follows;

x'=ϡ(y-x),y'=rx+-xz-y,z'=xy-bz,

where ϡ,r,b are constants. Perfom the stability analysis of this nonlinear system.

12. For the Holling-Tanner type Predator-Prey model

x'=x(1-x)-xa+xy,y'=y(ϒ-ϐyx),

where x(0)>0,y(0)>0 and a,ϒ,ϐ are positive constants. Find the equilibria and classify them.

13. Consider the Delayed Predator-Prey model

x'(t)=rx(t)-bx(t)x(t-Ϣ)-Ϗx(t)y(t),y'(t)=-cy(t)+ϐx(t)y(t),

with time delay Ϣ and positive constants r,b,Ϗ,c and ϐ. By choosing Ϣ as bifurcation parameter, check whether bifurcating periodic solutions occur around the equilibrium points or not.