## Finding a Hopf bifurcation from a cube equation

Let me first state what I found out this week. Forgive me if this result is trivial for you, but this result indeed helped me in my research.

*Lemma 1
Suppose we have a cube equation . The roots of this equation are a couple of purely imaginary roots and a single real root if and only if the following conditions are satisfied and .*

The proof is very easy and I am not going to prove this. Apparently, the conditions and are a necessary and sufficient condition for the cube equations to have a single real root and a couple of purely imaginary roots.

Okay, now let us move to the beginning why I bother with cube equations in the first place.

At the moment I am studying a mathematical model on HIV-1 virus with cell mediated immunity. I just started so I still don’t know all the representation of this model in the real life. But here is the model.

Yes, you are right, it is a five dimensional model. The variables are the state variables that are representing the density of uninfected cells, the infected cells, the density of virus, the CTL response cells and the CTL effector cells respectively. Why it is five dimensional, what are the parameters representing and other modeling questions, I still don’t know the answer. At the moment I am interested to find an equilibrium in this system that can undergo a Hopf bifurcation. An equilibrium can have a Hopf if the Jacobian of the system evaluated at this equilibrium has a purely imaginary eigenvalue. Actually this is the first time I handle a five dimensional system so it will be interesting.

Now we know why I need Lemma 1.

In my case, the situation is greatly simplified as all parameters are numbers except two which are and . From the paper that I read, I found that there are three equilibria in the system, so called . At the moment I am interested in the equilibrium because I suspected that this equilibrium can have a Hopf bifurcation and because the coordinate of this equilibrium has the following conditon . Let us have a look at the Jacobian matrix evaluated at .

As I said before, all the parameters are numbers except and . Also will be an expression that depends on and . As a result the matrix above will be depending on .

We are interested in what value of such that the matrix above has a purely imaginary eigenvalues. Even though the above matrix is not a block matrix but when we compute the eigenvalues, we can really separate the matrix above. We can obtain the eigenvalues of the above matrix by computing the eigenvalues of the following two matrices:

and

.

The eigenvalues of matrix is easy and it is real. So the problem really is in the computing the eigenvalue of matrix . Consider the matrix in a simpler notation.

.

The characteristic polynomial of the above matrix is

Therefore to find the value such that we have a Hopf bifurcation, we only need to solve the following conditions:

1. make sure that is positive and

2. solve .

I created this note (which is a part of our paper) so that I won’t forget what I had done. We don’t usually show this kind of computation but for me this little computation will be useful. Even though I uses software to compute and investigate the existence of Hopf bifurcation but it does not show the exact value of the parameter. Using algebraic approach I found the exact value.

Reference of the model

Yu, Huang, and Jiang, Dynamics of an HIV-1 infection model with cell mediated immunity, Communications in Nonlinear Science and Numerical Simulation, 2014, vol. 19.

Additional note: (Two dimensional case) Suppose the Jacobian matrix evaluated at the equilibrium is given below

.

Then the equilibrium has a Hopf bifurcation if the following conditions are satisfied:

1. and

2. .

## My first logistic bifurcation curve and a bit more

I don’t know since when, but I desperately want to learn how do I make a logistic bifurcation curve. And now, I finally able to create one. I do not want to forget this ability therefore I post this post, just in case I forget.

Logistic map is one of discrete dynamical systems, represented in the following:

.

It is basically a recursive sequence. It looks simple but it is highly not.When the sequence generated will always be inside the interval [0,4].More specifically, when , the sequence converges to a point. When is varied, does the solution still converge to some point?

In fact, when we vary past 3 we have the situation that is best explained through the following figure:

The figure is saying when is just past 3, the sequence is no longer converging to a point. This phenomena is called a period-doubling, as the sequence now tends to alternatingly approaching two points. However, this is the code to create the above figure in Matlab.

clear; clc;

a = 3:0.005:4; %interval of parameter

x0 = 0.5; %initial condition

N=500; %number of iterations

x(1) = x0; %the first entry is the initial condition

%computation of the orbit

figure(3);hold on;

for i=1:length(a)

for j=2:N

x(j) = a(i)*x(j-1)*(1-x(j-1));

end

y = x(301:end);

plot(a(i),y);

end

hold off

clear

Okay, the fitst task is done. Let us now learn how to make a movie in Matlab. I choose the topic ‘coweb’ as an example. Coweb arises in our topic, which is a discrete dynamical system.

Here is the video

The coweb created in the above movie uses . This exactly where the chaotic dynamic occurs. Therefore, we do not see any pattern there. The iteration that I used is only 90. However, we already see that the sequence makes such a complicated figure.

The code to create such a video is the following.

`clc;clear;`

myu = 3.9; %initial parameter

maxiter = 90; %maximum iteration

xo = 0.4; %initial condition

aviobj = avifile ( 'logistik2.avi', 'fps', 3 );%create a file logistik2.avi

x = 0:0.01:1; %range of x data

y1 = myu.*x.*(1-x); %range of y data to graph y=myu(x-x^2)

y2 = x; %range for y data to graph y=x

figure(1);hold on; %open figure window and hold on

axis([0 1.1 0 1.1]) % set the y and x axis

plot(x,y1); %plot the graph of y=myu(x-x^2)

plot(x,y2); %plot the graph of y=x

frame = getframe ( gca ); %command 1 to store what is inside figure 1

aviobj = addframe ( aviobj, frame ); %command 2

xbaru = myu.*xo.*(1-xo);

z2 = 0:0.01:xbaru;

z1 = xo;

plot(z1,z2);

frame = getframe ( gca );

aviobj = addframe ( aviobj, frame );

fprintf('the fixed point is %f',(myu-1)/myu);

for i=2:maxiter

if(xbaru > xo)

y1 = xo:0.01:xbaru;

else

y1 = xbaru:0.01:xo;

end

y2 = xbaru;

plot(y1,y2);

frame = getframe ( gca );

aviobj = addframe ( aviobj, frame );

xo = xbaru;

xbaru = myu.*xo.*(1-xo);

if(xbaru > xo)

y1 = xo:0.01:xbaru;

else

y1 = xbaru:0.01:xo;

end

y2 = xo;

plot(y2,y1);

frame = getframe ( gca );

aviobj = addframe ( aviobj, frame );

end

hold off;

H = (myu-1)/myu;

fprintf('please press any key \n');

aviobj = close ( aviobj );

%pause;

%close all;

However, I really really want to learn how to convert this avi video produced by Matlab to a gif image format. So if anyone know how do I do this please let me know.

## Your says