Open access peer-reviewed chapter

Bifurcation Analysis and Its Applications

By Canan Çelik Karaaslanlı

Submitted: November 17th 2011Reviewed: May 22nd 2012Published: September 19th 2012

DOI: 10.5772/50075

Downloaded: 6184

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,Ϗ),E1

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,Ϗ),E2
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)0E3

at the critical value Ϗ0where 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 fis a function mapping RnxR3Rn.A point x̣is called an equilibrium point if there is a specific ụϓRmsuch that

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

Suppose x̣is an equilibrium point (with the input ụ). Consider the initial condition x(0)=x̣, and applying the input u(t)=ụ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, rand Kare positive constants, Kis 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=0and 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 xto 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 ϓ0is a constant. Hence, the solution is growing if f'(x̣)>0and decaying if f'(x̣)<0.As a result, the equilibrium point is stable if f'(x̣)<0,unstable if f'(x̣)>0as 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 Iwhere the equilibrium point x̣ϓI. Then the equilibrium point x̣is locally stable if f'(x̣)<0and 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

fṆ=0,whichyieldstwopointsṆ=0,K.

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

f'(N)=r(1-2NK)atṆ=0,f'(0)=r>0,soitisunstableatṆ=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̣,ỵ)is a steady state (equilibrium point), i.e.,

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

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

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

where uand vare understood to be small as u1and v1.Is is natural to ask whether uand vare growing or decaying so that xand ywill 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,ỵ+v)=fx̣,ỵ+fxx̣,ỵu+fyx̣,ỵv+higherorderterms...=fxx̣,ỵu+fyx̣,ỵv+higherorderterms....

Similarly,

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

Since uand vare 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 uand v,

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

where the matrix fxfygxgyis called Jacobian matrix Jof the nonlinear system.The above linear system for uand vhas 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̣,ỵ)of the differential equation is stable if all the eigenvalues of J, the Jacobian evaluated at (x̣,ỵ)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 Jare negative, or have negative real parts.

Unstable: A critical point is unstable if at least one eigenvalue of the jacobian matrix Jis 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̣,ỵ)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)],E4
y(t)=x(t)[2-y(t)].

The equilibria are the points (x̣,ỵ)=(0,0)and (x̣,ỵ)=(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-120which 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-1thus the point is locally unstable (as

ϙ1=1andϙ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,ϙ),E5

where f:Rn+1Rnis 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 ϙ0is a bifurcation point of equation (5).

3.1. Bifurcation in one dimension

We may assume that n=1so 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 ϙ0is a bifurcation point if ϙ1(ϙ0)=0.Hence, bifurcation points (x0(ϙ),ϙ)are solutions of

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

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

3.1.1. Saddle-node bifurcation

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,foraisreal.E6

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|0as t(linear stability); for x=-a,|x|0as t(linear stability).

As sketched in the “ bifurcation diagram” below, therefore, the saddle node bifurcation at a=0corresponds 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,forx,a,breal.

Again, aand bare 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=0corresponds 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,fora,breal.

As usual, aand bare external control parameters. Steady states, for which x'=0are 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̣2and x̣3only exist when a>0if b>0and for a<0if 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̣2and 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,Ϣ,E7
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 Ϣ*=0gives purely imaginary eigenvalues. System (7) is rewritten as follows;

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

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

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

is the Jacobian matrix evaluated at origin.

Theorem (Hopf bifurcation theorem)

Let f1and g1, in system (8) have continuous third order partial derivatives in xand 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,ϐ00such that the eigenvalues cross the imaginary axis with nonzero speed, i.e., dϏdϢ|Ϣ=00.

Then in any open set Ucontaining the origin in R2and for any Ϣ0>0,there exists a value Ϣ̣,|Ϣ̣|<Ϣ0such 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-Ϗ,E9
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-Ϗ-1and 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 Ϗ=0is the bifurcation point for this differential equation system.

Example: Consider the two dimensional system

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

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

J=Ϗ1-1Ϗ

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

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

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

4. Center Manifold teorem

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

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,n0number of eigenvalues with Reϙ=0and n-number of eigenvalues with Reϙ<0.Let Tcbe the eigenspace on imaginary axis corresponding to n0eigenvalues. The eigenvalues on the imaginary axis Reϙ=0are called the critical eigenvalues as on the eigenspace Tc.And suppose the function Ϥtdenote 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 Ccenter manifold Wlocc0such 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 Asare are the blocks in the canonical form whose diagonals contain the eigenvalues with Reϙ=0and Reϙ<0; respectively) restricted to the center manifold are given by

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

And the manifold Wloccis 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,E12
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=0and the center manifold reduction takes the form

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

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

Example: Consider the two dimensional system

x'=x2y-x5,E13
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=1and b=0and

x'=x4+O(x5).

For this reduced equation, x=0is 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ϓR3E14

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ϧ1is 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),E15

where x(t)=(x1(t),x2(t),x3(t))TϓR3,and LϚ:CR3,f:R×CR3are 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Ϝ,E16

where ϕ(ϖ)=ϕ(ϖ,0).Then A(0)and A*are adjoint operators. Suppose that q(ϖ)and q*(s)are eigenvectors of Aand A*corresponding to iϧϢkand -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 C0at Ϛ=0.Define

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

On the center manifold C0, we have

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

where zand ẓare local coordinates for center manifold C0in the direction of qand q̣*.Note that Wis real if xtis real. We consider only real solutions. For the solution xtϓC0,since Ϛ=0and (15), we have

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

where

g(z,ẓ)=q̣*(0)F0(z,ẓ)=g20z22+g11zẓ+g02ẓ22+g21z2ẓ2+....E18

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

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

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

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

and we evaluate g(z,ẓ).

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

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

where

H(z.ẓ,ϖ)=H20(ϖ)z22+H11(ϖ)zẓ+H02(ϖ)ẓ22+....E20

Note that on the center manifold C0near to the origin,

W̥=Wzz̥+Wẓẓ.E21

Thus we obtain,

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

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

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

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

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

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))ϓR3is a constant vector. From the definition of Aand (22), we obtain

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

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'sand 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)}ϧϢkE27

and we state this as in the following theorem.

Theorem: Ϛ2determines 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. ϐ2determines the stability of the bifurcating periodic solutions; bifurcating periodic solutions are stable if ϐ2<0,unstable if ϐ2>0. T2determines 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).E28

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),E29

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=ϧ12and g'(z1)=1.1944×10-4>0which 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.5663and also Hopf bifurcation occurs at Ϣ=Ϣ0=1.5663as 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, ϐ2and T2as

Ϛ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.5663is 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<Ϣ0and 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>Ϣ0sufficiently close to Ϣ0which 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 y0be an equilibrium point of the equation

dydt=f(y).Let g(y)be a positive function. Deduce the stability of y0as 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 ras 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*=0is 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 rthere are additional fixed points. For which values of rdo 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 rand 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. r2=x2+y2.r'=0
  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,bare 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)>0and 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,Ϗ,cand ϐ.By choosing Ϣas bifurcation parameter, check whether bifurcating periodic solutions occur around the equilibrium points or not.

© 2012 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Canan Çelik Karaaslanlı (September 19th 2012). Bifurcation Analysis and Its Applications, Numerical Simulation - From Theory to Industry, Mykhaylo Andriychuk, IntechOpen, DOI: 10.5772/50075. Available from:

chapter statistics

6184total chapter downloads

1Crossref citations

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Model-Based Adaptive Tracking Control of Linear Time-Varying System with Uncertainties

By DongBin Lee and C. Nataraj

Related Book

First chapter

Modeling the Physical Phenomena Involved by Laser Beam – Substance Interaction

By Marian Pearsica, Stefan Nedelcu, Cristian-George Constantinescu, Constantin Strimbu, Marius Benta and Catalin Mihai

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More about us