Archive

Posts Tagged ‘bifurcation’

Finding a Hopf bifurcation from a cube equation

September 27, 2014 1 comment

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 x^3 + ax^2 + bx + c=0. 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 b>0 and c=ab.

The proof is very easy and I am not going to prove this. Apparently, the conditions b>0 and c=ab 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.

\dot{x} = \lambda - dx - \beta x v
\dot{y} = \beta x v - ay - p y z
\dot{v} = k y- u v
\dot{w} = c y w - c q y w - b w
\dot{z} = c q y w - h z

Yes, you are right, it is a five dimensional model. The variables x, y, v, w, z 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 d and b. From the paper that I read, I found that there are three equilibria in the system, so called E_0, E_1, E_2. At the moment I am interested in the equilibrium E_1 because I suspected that this equilibrium can have a Hopf bifurcation and because the coordinate of this equilibrium has the following conditon w=z=0. Let us have a look at the Jacobian matrix evaluated at E_1=(x_1,y_1,v_1,0,0).

J = \left( \begin{array}{ccccc} -d-\beta v_1 & 0 & -\beta x_1 & 0 & 0\\\beta v_1 & -a & \beta x_1 & 0 & -p y_1 \\ 0 & k & -u & 0 & 0 \\ 0 & 0 & 0 & cy_1-c q y_1 - b & 0 \\ 0 & 0 & 0 & c q y_1 & -h \end{array} \right)

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

We are interested in what value of d,b 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:

A = \left( \begin{array}{ccccc} -d-\beta v_1 & 0 & -\beta x_1 \\ \beta v_1 & -a & \beta x_1 \\ 0 & k & -u \\ \end{array} \right) and

B = \left( \begin{array}{ccccc} - cy_1-c q y_1 - b & 0 \\ c q y_1 & -h \end{array} \right).

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

A = \left( \begin{array}{ccccc} a & b & c \\ d & e & f \\ g & h & i \\ \end{array} \right).

The characteristic polynomial of the above matrix is

p(\lambda) = \lambda^3 - (a+d+g) \lambda^2 + (bd-ae+fh-ei+cg-ai) \lambda + \det(A)

Therefore to find the value d,b such that we have a Hopf bifurcation, we only need to solve the following conditions:
1. make sure that (bd-ae+fh-ei+cg-ai) is positive and
2. solve - (a+d+g) \times (bd-ae+fh-ei+cg-ai) = \det(A).

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

B = \left( \begin{array}{ccccc} p & q \\ r & s \end{array} \right).

Then the equilibrium has a Hopf bifurcation if the following conditions are satisfied:
1. \det(B) > 0 and
2. trace(B) = 0.

 

 

My first logistic bifurcation curve and a bit more

May 27, 2011 5 comments

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:

x_{n+1} = \mu x_n (1-x_n).

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

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

The figure is saying when \mu 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 \mu=3.9. 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.