Programming and writing about it.

echo $RANDOM

Tag: math

Doing Math with Python: Chapters 5 and 6 in early access

I am excited to share that the fifth and sixth chapters are available as part of the early access of my book Doing Math with Python.


Chapter 5: Sets and Probability

This chapter starts off with how to create a set and demonstrating the common set operations. Utility of the different set operations are demonstrated via simple applications. For example, Cartesian product is used to write a program to simulate an experiment to calculate the time period of a simple pendulum of different lengths and at places with varying gravity. Union and intersection operations are applied to finding the probability of events.

The chapter then moves onto discussing how to generate uniform and non uniform random numbers, and using them to simulate scenarios such as a die roll and a fictional ATM which dispenses dollar bills of different denominations with varying probability.

One of the challenges at the end discusses drawing venn diagrams.

Chapter 6: Drawing shapes and Fractals

This chapter is logically divided into two parts. The first part introduces the reader to matplotlib patches which allows drawing geometric shapes (circles and polygons), followed by matplotlib’s animation API which allows drawing animated figures. The trajectory of a projectile motion discussed elsewhere in various contexts is animated combining both these things.

The second part of the book introduces the concept of geometric transformation. Combining that with the knowledge of generating random numbers learned earlier in Chapter 5, the reader will learn how to draw fractals such as the Barnsley Fern.

The challenges at the end gives the opportunity for the reader to explore the Sierpinski triangle and Henon’s function.

Trying out the programs

Using the Anaconda distribution (Python 3) should be the easiest way to try out all the programs in the book. You will need matplotlib, sympy 0.7.6 and matplotlib_venn to try out the programs. An installation guide will be available online soon.

Stay updated

I am working on the last chapter for the book. You can stay updated on the book via various channels:

Blog posts:

Facebook page:

G+ Community:


If you are interested in taking a look at a sample copy, I can try to get a sample for you to look at the current pre-released version of the book. Please feel free to get in touch.


Doing Math with Python: Stay Updated

I am reaching the final stages of my new book. Here are few ways to stay updated about the book:


Blog posts:

Facebook page:

G+ Community:


If you are an educator/teacher, I can also try to get a sample for you to look at the current pre-released version of the book.

Doing Math with Python: Two more chapters in Early Access

I am excited to share that the third and fourth chapters are available as part of the early access of my book Doing Math with Python.


Chapter 3: Describing Data with Statistics

As the title suggests, this chapter is all about the statistical measures one would first learn in high school – mean, median, mode, frequency table, range, variance, standard deviation and linear correlation are discussed.

Chapter 4: Algebra and Symbolic Math with SymPy

The first three chapters are all about number crunching. The fourth chapter introduces the reader to the basics of manipulating symbolic expressions using SymPy. Factorizing algebraic expressions, solving equations, plotting from symbolic expressions are some of the topics discussed in this chapter.

Trying out the programs

Using the Anaconda distribution (Python 3) should be the easiest way to try out all the programs in the book.

Announcing early access for “Doing Math with Python”

Early Access promotion:

No Starch Press is running a promotion tomorrow (29/01/2015 – PST)  that will offer 40% off all Early Access titles on, including Doing Math with Python. The discount code is BRIGHTANDEARLY. Currently two chapters are available and more should be up in the coming days!

I am excited to write that my new book “Doing Math with Python” to be published by No Starch Press is now available via their “Early Access” program – which means if you now buy the book, you will get the chapters as and when they are available and also, the chapters may need more polishing.


The book uses Python 3 exclusively and the Appendix A covers setup and installation instructions for Python 3 and the libraries used in the book. However, that is not yet available. Hence, at this stage, the easiest way would be to use a distribution such as Anaconda on Windows, Linux or Mac OSX (untested, I don’t have access to the OS).

Circles and Infinity

As we go outword from the inner circle, the gaps between the lines joining the centre to the circumference widens:

This essentially means that even though you can theoretically you can have an infinite number of lines on the inner circle, they are never enough to cover the outer circle(s). I was made to ponder on this fact while watching a BBC documentary titled “Dangerous Minds”.

Here is the MATLAB code (works in GNU Octave):

clear all;close all

% radius
for r=1:10


    for theta=0:360


    % plot the circle
    plot(cx,cy,'r.'); hold on

    % lines matrix
    lines=[lines cx cy];

    % plot the lines
    for i=1:size(lines,1)
        plot([lines(i,1) lines(i,3)],[lines(i,2) lines(i,4)],'b-');
        hold on


And here is another thing for you to think about, also mentioned in the documentary: A circle can be thought of as having infinite points on its circumference. The bigger circle looks like having more points on the circumference, right?

Now consider these two circles are concentric: so for every point on the outer circle, a line joining it to the centre intersects the inner circle at a point.  So, is there really more points on the circumference of the bigger circle? Time to blame infinity.

Eigen values, vectors and Fibonacci sequence

I came across application of using Eigen values and Eigen vectors to find the N’th Fibonacci number very easily in Prof. Gilbert Strang’s lecture series.

Fibonacci numbers

Fibonacci numbers are a sequence of numbers such that every term is a sum of the previous two terms: F_{n+2} = F_{n} + F_{n+1}. The first few terms of this series are 0,1,1,2,3,5,8,13..

Linear Dynamical System

Consider a N-dimensional vector \vec{X}. Assume that the variation of \vec{X} with respect to time, when it varies in discrete steps may be represented as: X_{t+1}=AX_{t}, where A is constant matrix  Such a system is called a Linear Dynamical System. For example, X_{n+1}=2X_n is a simple Linear Dynamical System.

Fibonacci Number sequence as a Linear Dynamical system

Let us now make an attempt to model the Fibonacci number sequence as a Linear Dynamical System. This is our first step towards this beautiful application of Eigen values and vectors.

Consider the following equations:

F_{n+2}=F_n + F_{n+1}                                                                                                                                       F_{n+1}=F_{n+1}

The second equation is written to enable us to write the above as a Linear dynamical system. Now, can we write the above equations as follows:

which is of the form U_{n+1}=AU_{n}, where A is the constant matrix as above.

The matrices S and \Lambda

The matrix S is said to be the Eigen vector matrix, i.e S= [x_1 x_2 .. x_n], considering x_1, x_2, , x_n are the n Eigen vectors of the constant matrix A, defined above. The \Lambda matrix is a diagonal matrix defined as follows:

where \lambda_{1}, \lambda_{2}, \lambda_{3}, .. \lambda_{n} are the Eigen vectors of the same matrix.

From the basic equation of Eigen values and vectors Ax= \lambda x, we can write the above as AS=S\Lambda. From this equation, we can write :                                                                                                                                                                                                  A = S \Lambda S^{-1} , AA = (S\Lambda S^{-1}) (S\Lambda S^{-1}), AA = S\Lambda (S^{-1}S) \Lambda S^{-1} = S\Lambda ^{2}S^{-1}, A^2 = S\Lambda ^{2}S^{-1} , A^n = S\Lambda ^{n}S^{-1}

This is an useful and important relationship. How? Earlier, we wrote:   U_{n+1} = AU_{n}, U_{1} = AU_{0} , U_{2} =AU_{1} = A(AU_{0}) = A^{2}U_0 U_n = A^{n}U_{0}

Now, substitute A^n  by S\Lambda ^{n}S^{-1}U_0 can be written in a more simplified manner, by writing it as a combination of its Eigen vectors, U_0=Sc where c is a coefficient matrix, from which we can write c=S^{-1}U_0. Thus we have:

U_n = S\Lambda ^{n}S^{-1}Sc

U_n = S\Lambda ^{n}c

So what does the last equation say? It says, the current state (n^{th} state) of our model is expressed as matrix obtained by the product of the Eigen vector matrix S, the Eigen value matrix \Lambda, raised to the power k and the matrix c.

Let us now make use of this equation to get the Nth Fibonacci number.

Eigen Values and Vectors of A

The Eigen values of A can be shown to be: \lambda_{1}= \frac{1+\sqrt{5}}{2} and \lambda_{2}= \frac{1-\sqrt{5}}{2}. The Eigen vectors of A turns out to be x_1 = \frac{1+\sqrt(5)}{2} and x_2 = \frac{1-\sqrt(5)}{2}.

S and \Lambda

The coefficient matric c
As derived earlier,  c=S^{-1} u_0. We already have S, S^{-1} turns out to be:

The first two Fibonacci numbers are 1,1, so we can take u_0 =\begin{bmatrix} 1 \\ 1 \end{bmatrix}. Hence:

Finding the Nth Fibonacci Number

Now, let’s rewrite the earlier equation U_n = S\Lambda ^{n}c, by substituting the values such that it matches our Fibonacci model:

We can rewrite the above equation as follows:

So, we have it! Just substitute the value of ‘n’ and we get the F_{n+1} and F_{n+2} Fibonacci numbers. For example, let n=0 and see the values of F_1 and F_2, they should come out to be 1, 1 respectively. This piece of Python code prints the first ‘n’ Fibonacci numbers using the above formulae:

# Print fibonacci numbers by forming it as a linear dynamical
# system
# Amit Saha (

import math
import numpy as np

# Eigen vectors of A
eig_vec1= (1.0 + math.sqrt(5))/2.0
eig_vec2= (1.0 - math.sqrt(5))/2.0

# coefficients
c_1 = ((5.0 + math.sqrt(5.0))/10.0)
c_2 = ((5.0 - math.sqrt(5.0))/10.0)

# matrices
mat1 = np.matrix([[eig_vec1],[1]])
mat2 = np.matrix([[eig_vec2],[1]])

# set n accordingly to get the first (n+1) fibonacci numbers
for n in range(0,n):
    sum_mat = (eig_vec1**n)*c_1*mat1 + (eig_vec2**n)*c_2*mat2
    print float(sum_mat[1])

Parting notes

A very interesting thing to note here is that the second term of the last equation, contains the term \frac{1-\sqrt{5}}{2}. The absolute value of this term t is between 0 and 1, hence as ‘n’ increases its value goes closer and closer to 0 and almost vanishes. Thus, with large values of ‘n’, the value of the above expression can be assumed to be determined by the first term. This has an interesting implication, for which I shall direct you to this awesome post.

Also check Project Euler problem #2: