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.




Math on the web 2

September 11, 2014 Leave a comment

This is a second edition of Math on the web, in which I try to list mathematical articles around the internet. I found these articles from googling, from social media networks and many more sources. I will write website links that relate to mathematics and its applications that I found from May 2014 until now. I hope you find them interesting because these are really very subjective.

  1. @TopologyFact : Notes on hyperbolic geometry
    The above link will bring you to a kind of ebook about hyperbolic geometry. To be honest, I have not read it. I just downloaded it and it now sits in my computer hard drive. I got this because in case I want to learn about hyperbolic geometry.
  2. is a webpage that (I think) is similar to the (late) google reader. It can bring you the latest information you are interested in and even it can bring you information from social media like Facebook, Google Plus and Twitter. The link above is a paper created by Nick Chater that brings us information about Mathematical Education. I believe the information will be very useful for mathematics educators like me and you.
  3. @approx_normal: A few awesome bumper stickers here
    The link will bring you to various bumper stickers that are very cool when you are mathematicians. I would buy it if there is someone selling these, especially with the Lorrentz Attractor one.
  4. @kdnuggets: Numeric matrix manipulation: cheat sheet for MATLAB, Python NumPy, R, and Julia #DataScience I am a fan of mathematical softwares. I could sit in my office for hours trying to understand new language of mathematical softwares. This link will bring us the cheat sheet of some most popular languanges. It is not a complete cheat sheet though as it is only discussing the basic operation of matrices., which can be very useful for people who just begin to learn the language.
  5. @pickover Teachers, introduce your students to the wonders of Hilbert Curves at this awesome website.
    @pickover is a twitter acount of a mathematician named Cliff Pickover . He is a well known author who has published more than 50 books. Cited from his website, his books explore topics ranging from computers and creativity to art, mathematics, parallel universes, Einstein, time travel, alien life, religion, dimethyltryptamine elves, and the nature of human genius.
    The link above will bring a very introductive article about Hilbert Curve.
Categories: Teachings

Math on the web 1

May 23, 2014 Leave a comment

The idea just came out of the blue. I like mathematics and I am using twitter. I am amazed that there are so many mathematicians using twitter. I am soo fortunate that there are so many mathematics articles around the web that are very useful, interesting and opening our minds. So, why not save and keep those interesting articles (according to me) .  These articles are mathematics articles that I read in the past few months.

1. @MrHonner: Evolutionary game theory in the classroom by @drvinceknight #mathchat #math

@MrHonner is a twitter account owned by a mathematician Patrick Honner. I followed him and at some stage his tweet opened my mind to learn something new in mathematical teaching

2. @MrHonner: “It’s not how big your class is, it’s what you do with it” Interesting take on class size #edchat

Again, this was a tweet from @MrHonner that brought me an article about the class size vs the quality of teaching. Since we are all teacher, I am sure that this article will be of interest.

3. @ISACalculus: Daylight Savings is mathematically illogical. Dilation not translation. Thanks @MrHonner for clearing this up.

This is a light article yet interesting about Daylight Saving. I am now living in Indonesia, a very nice tropical region which does not have Daylight Saving. I once lived in Australia where Daylight Saving occurs every November – March (the other way around from Europe). I find it Daylight Saving is very interesting and I like this article very much.

4. @standupmaths: You need some maths inequalities? This guy has ≥ you’ll ever need:

Very useful article for mathematics students who wants to compete in mathematics olympiad. But let me give a credit to @standupmaths first. This account belongs to Matt Parker who is (of course) a mathematician. And he is also a stand-up comedian. If I get the chance and the universe allows me, I would invite him to my university. Anyway, the title of the article explains a lot. It would give us a summary of inequalities you have ever (or never) heard. So, grab it before the link is dead.

5. I am so sorry that I do not put who is the one tweeted this link cause I have forgotten. However, the link will bring you to an article about fractal. It is a very nice introduction to fractal and you’ll love it the way the author presented the article. Not only did he provide figures and illustrations but also he provided animations and simulations on how to make a fractal shape. 

6. Again I have forgotten who was the one who tweeted/retweeted this link. This is my last article about mathematics on the web. I think it is quite pretty to end this post with this article cause this article is about how math can give you the best job (in terms of salary). The article said in the beginning: “Another day, another reason to get better at math,” which for me, is true.

As a conclusion, I’ll let you read the articles. I am pretty sure I am going to write math on the web 2 next month or after that.

Categories: Teachings

Random Numbers in Matlab

May 23, 2014 3 comments

In the past two years, I have been supervising bachelor degree students for their final projects. Unfortunately, most of the projects were related to time series, forecasting, stochastics process, financial mathematics and many topics related to application to statistics in finance. To be honest, this is not my strong suit. As a result, we struggled when reading and using many statistical techniques especially the one in time series. With this notes, I hope to understand and to keep remember those techniques. Moreover, as I am teaching Time Series Analysis this semester, I believe this post will be very useful for my students.

In this post, I concentrate on how to generate random numbers with certain distributions in Matlab and its application in my time series lecture. I know that we can easily google on how to do these, but as always, it would be very nice to have this ready in my blog. When I forget to do this, I don’t have to google it again and save myself a few minutes.

Generating random numbers and introduction to subplot, histogram, state commands in matlab

To generate random numbers in Matlab with certain distributions, I  need to type

x1 = normrnd(mean,std,[m,n]); %normal distribution
x2 = binornd(N,p,[m,n]); %binomial distribution
x3 = exprnd(lambda,[m,n]); %exponential distribution
x4 = gamrnd(a,b,[m,n]); %Gamma distribution
x5 = rand(m,n); %uniform distribution
x6 = poissrnd(lambda,[m,n]); %Poisson distribution

The commands above returns an m x n matrix containing pseudorandom values drawn from normal, binomial, exponential, gamma, uniform and Poisson distributions respectively. Some usefull commands related to random number generator in Matlab are

subplot(2,2,1) %top left figure
hist(x1) %histogram of normal random numbers
subplot(2,2,2) %top right figure
hist(x2) %histogram of binomial random numbers
subplot(2,2,3) %bottom left figure
hist(x3) %histogram of exponential random numbers
subplot(2,2,4) %bottom right figure
hist(x4) %histogram of Gamma random numbers

The above commands return histograms of random numbers generated by the previous commands.  The output of those commands is the following figure.


Another useful command regarding random number generators is the ‘state’ option. This is very usefull if we want to repeat our computation that involves random number. The following code will help us understand how to use the state command.

theta = -0.8;
phi = 0.9;
mu = 0;
sigma = 0.1;
%-----generating random numbers-------------------------
e = normrnd(mu,sigma,[100 1]);
%-----generating MA(1) and AR(1) process----------------
y(1) = 0;
z(1) = 1;
for i = 2:length(e)
y(i) = e(i) - theta*e(i-1);
z(i) = phi*z(i-1) + e(i);
%-----plotting MA(1) and AR(1) process----------------

Using the above command, especially when we write the ‘ state’ option in line 8, it allows us to repeat our computation later. We can show the result of our computation tomorrow exactly the same as our computation we conduct today, even though our code involves a random number.  The result of the above command is plotted in the following figure.


Plotting the auto correlation function and introduction to bar command

In time series analysis, when we have a time series, it is common to plot the sample auto correlation function (ACF) to be compared to time series models we have in the textbooks. The thing is today, our students now rely heavily on statistical softwares to plot the ACF so that they forget how to plot it in the beginning.  Recall that the sample ACF of the observed series Y_1, Y_2, \dots, Y_n is defined as

(1)——–\displaystyle{r_k = \frac{\sum_{t=k+1}^{n}(Y_t - \bar{Y})(Y_{t-k}-\bar{Y})}{\sum_{t=1}^{n}(Y_t - \bar{Y})^2}}

The following code shows a MATLAB code that will plot the ACF of both time series generated by the code above.

sumY = 0;
sumZ = 0;
for i = 1:length(y)
sumY = sumY + y(i);
sumZ = sumZ + z(i);
ybar = sumY/length(y);
zbar = sumZ/length(z);
sum2Y = 0;
sum2Z = 0;
for i = 1:length(y)
sum2Y = sum2Y + (y(i)-ybar)^2;
sum2Z = sum2Z + (z(i)-zbar)^2;
for k = 1:length(y)
sum3Y = 0;
sum3Z = 0;
for t = k+1 : length(y)
sum3Y = sum3Y + (y(t)-ybar)*(y(t-k)-ybar);
sum3Z = sum3Z + (z(t)-zbar)*(z(t-k)-zbar);
ry(k) = sum3Y/sum2Y;
rz(k) = sum3Z/sum2Z;

The last two commands above will plot the ACF of time series y and z in the bar style as it is always for the ACF, as shown in the following figure.




It doesn’t show, really, the supposed to be acf of MA(1) (above) and AR(1) (below) process. Perhaps this is I have not learned  the sample acf measurement (1).

How to export/import data between Matlab and Excell

Finally, we are going to write a code in Matlab on how to export or import data to/from Excell in Matlab. Again, this is to save myself a few minutes rather than googling it around the internet.

%--------------write data to excell-----------------------
y = transpose(y);
z = transpose(z);
filename = 'ma_ar_data.xlsx';
%-------------read data from excell----------------------
filename = 'ma_ar_data.xlsx';
sheet = 1;
xlRange = 'A1:A100';
dataY = xlsread(filename, sheet, xlRange);

The transpose command above puts the variable y generated before to be a column vector.


In this post, I have created various codes so that I am able to remember

  • generating random numbers
  • using subplot
  • using histogram and bar
  • exporting/importing data to/from excell
  • plotting the ACF


Liapunov exponent computation

October 4, 2013 Leave a comment

Liapunov exponent method is one of the clever ways to investigate a certain dynamical system to be chaotic. In this post, I shall try to compute Liapunov exponent of one famous discrete dynamical system, i.e. the logistic map. After this, I shall try to compute the Liapunov exponent of a continuous dynamical system.

Actually, I write this post because I want to know how to numerically compute the Liapunov exponent. I often read a book or a website such as wolfram.alpha that shows the Liapunov exponent but it didn’t explain how to get this numerically. After I am able to post this post, and when I forget how to compute this exponent, I’d know where to look.

Suppose we have an initial condition x_0 and a nearby initial point x_0 + \delta_0, where \delta_0 is an initial difference between initial points, which by the way will be extremely small. These two nearby points will be iterated by some map. Let \delta_n be the difference between these two points after n iterations. Consider the following equation

(1) …  |\delta_n| =|\delta_0 |e^{n \lambda}.

The coefficient \lambda is called the Liapunov exponent. We want to compute this coefficient. Note that \delta_n = f^n(x_0)-f^n(x_0+\delta_0), where f is the given map. Taking logarithm of (1), we get

(2) … \lambda = \frac{1}{n} \ln \left|\frac{\delta_n}{\delta_0}\right|=\frac{1}{n} \ln \left|\frac{f^n(x_0)-f^n(x_0+\delta_0)}{\delta_0}\right|=\frac{1}{n} \ln \left|{(f^n)'(x_0)}\right|

as \delta_0 \to 0. The term inside the logarithm can be simplified by using the chain rule: (f^n)'(x_0) = f'(x_{n-1}) f'(x_{n-2}) \dots f'(x_0). Thus (2) becomes

(3) … \lambda = \frac{1}{n} \sum_{i=0}^{n-1} \ln |f'(x_i)|.

When n \rightarrow \infty, \lambda is defined as the Liapunov exponent. When this exponent is positive, it is a sign of a chaotic attractor.

Let us consider the logistic map that was in my previous post here. Let us fix the parameter value r. Starting at a random initial point, we will iterate the initial point until it reaches the steady state (which will be a fixed point, periodic or aperiodic), say 300 iterates. Now we compute a large number iterates, say 10000 and compute (3) starting from the 301st point. Repeat this procedure for the next r. Then we plot the Liapunov exponent with respect to r. We shall find a figure as follows.


As we can see from the above figure. The map is chaotic in some values of parameter between 3.5 to 4, in which the Liapunov exponent is positive.

The code to produce the above figure is the following

tic; %to compute how long is the computation
clear;clc; %clear screem and clear all the previous variables
a = 3:0.001:4; %parameter of the logistic map
N = 300; %number of iterations until it is stable
N1 = 10000; %number of iterations to compute the Liapunov exponent
y = zeros(length(a),1); %declares variable for Liapunov exponents
for i=1:length(a)
x(1) = rand; %random initial conditions%
for j = 2:N
x(j) = a(i)*x(j-1)*(1-x(j-1));
sum = 0;
for j = 1:N1
x(N+j) = a(i)*x(N+j-1)*(1-x(N+j-1));
sum = sum + log(abs(a(i)-2*a(i)*x(N+j)));
y(i) = sum/10000;
fprintf('%d %d\n',i,length(a));
plot(a,0); hold on;
plot(a,y); %plot Liapunov exponent against parameter
toc; %display how long computation is

A bit of mathematics out of origami

February 27, 2013 1 comment

Dodecahedron Calendar

Last week on Saturday, we went to SMAK Hendrikus, one of senior High-schools in Surabaya. We were there to give a seminar for high-school teachers. The theme of the seminar is “the Amazing Math.”  Personally, I always think that Mathematics is amazing and I really want to convince this to other people.

Apart from the seminar, we were planning to give some souvenir for the participants. After consulting with my colleagues, we decided to give them a rhombic dodecahedron calendar. If there is anyone of you who has never heard of this kind of calendar, you can see it in the picture below, or you can follow this link to download the 2013 calendar along with the instruction on how to fold this model.rhombic_calendarWe were excited to know that this was actually a great gift for them. After reading the instruction we then realized the paper needed to make this calendar has to have a ratio which is normally a ratio of an A4 paper. Unfortunately, we already bought many origami papers which are all square.


The Lichtenberg Ratio

The Lichtenberg ratio is the ratio used in all of A-series paper such that A0, A1, et cetera. This ratio was found by Georg Christoph Lichtenberg. In the A-series paper the ratio between the length and the width of the paper is always \sqrt{2}:1. The unique property of this ratio is the following fact. Suppose the original paper is cut into two equal pieces (smaller pieces) such that the width of the original paper becomes the length of the smaller paper and the width of the smaller paper is half of the length of the original paper. Then the ratio between the length and the width of the original paper is the same as that of the smaller paper.the ratio between the length and the width stays the same.


We go back to our problem above. We do not have such a paper, but we have a bunch of papers whose dimension are all square. We then notice the following theorem that actually solves our problem.

Theorem (Construction a Lichtenberg ratio paper from a square)
Suppose we have a paper whose dimension is a \times a. If we follow the instructions below then we have a paper whose dimension follows the Lichtenberg ratio.


1. Draw a diagonal line as in the above figure.


2. Fold the paper along the diagonal line as the above right figure.


3. Mark the meeting point between the bottom right corner and the diagonal line as O in the above figure. then draw a vertical line. Cut the paper along this line.


4. We have a paper whose dimension is the Lichtenberg ratio as in the above figure.

Proof of the theorem

Consider the following figuresquare4

Let the square be the square ABCD as in the above figure. Assume the square has been folded and unfolded according the above instruction such that E, F, O and G are well defined. Suppose x is equal to AE, c is equal to CO and b is equal to CF. We are going to prove that

a : x = \sqrt{2} : 1.

First notice that AO is equal to a. Using the triangle ABC and the Pythagoras theorem, we get c = a\sqrt{2}-a.

Using the triangle OFC and the fact that the angle OCF is 45 degrees, we get b = a - a/\sqrt{2}.

Using the fact that the triangle OGC is an isosceles triangle we then obtain x=a-2b=a/\sqrt{2}.

This proves the theorem.

Categories: math puzzle, Teachings

Files output in C++

September 3, 2012 Leave a comment

All the programs that we normally use inputs from keyboard and give outputs on the screen. However, in a large scale computation, this situation is not effective. Thus, I have learnt how to produce outputs to a file, so that it can handle a large amount of data. I learned this back in 2007 and I have already forgotten about this code. Now, if I forget about this I know where to look.

Below I wrote a bit of source code that has to be compiled in C++. This program evaluates the exponential of a number from 0 to 3 using four different ways. The first three are using the Taylor polynomial of degree 1, 2 and 3 and the fourth way is using the built in exponential function from C++ itself. The output of this program is a file called “data.dat” that consists of five columns of numbers that can be used for plotting or other programs.

#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;
double taylor(double x, int n);
int factorial(int n);
int main()
double number;
ofstream out_stream;"data.dat");
for (number = 0 ; number <= 3.01 ; number = number + 0.01) {
out_stream << number << " ";
out_stream << taylor(number,1) << "\t";
out_stream << taylor(number,2) << "\t";
out_stream << taylor(number,3) << "\t";
out_stream << exp(number) << "\n";
return 0;
double taylor(double x, int n)
double y=1;
for(int k=1; k<=n; k++)
y = y + pow(x,k)/factorial(k);
return y;
int factorial(int n)
if (n <= 1)
return 1;
return n*factorial(n-1);
I then used gnuplot to plot the output of the program. The plot is the following
where column 1 is the x-axis, while columns 2,3,4 and 5 are the evaluation of exponential of x-axis using four different ways. We can see that Taylor polynomials converge to exp(x).