A Candle in the Dark

A look on science, politics, religion and events

Archive for the ‘mathematics’ Category

Abstruse Goose Puzzle

with 9 comments

You might have come across Abstruse Goose before. It’s simply a wonderful webcomic, along the lines of xkcd, but much more technical. Each comic strip usually contains an allusion to some physical or mathematical concept and ends with a humorous punchline.

Today’s comic had a simple puzzle. And, as I was curious about where it was leading, and had time to kill, I took a shot at it.

Clue 1

And how does one go solving this? It is fairly easy. In fact, all you need is to know is the first n digits of e and after that, it’s simply string manipulation and prime number checking.

I’ve embedded the C++ code which I used to solve the puzzle at the end of the post, below the fold.

It’s easier to find the first 100,000 digits of e using a quick search in google than to generate it. Once you have the string, you just need to check if every set of 10 consecutive digits is prime. Turns out, that the first 10 digit prime number found in consecutive digits of e is 7427466391.

That’s Clue #1 done. Proceeding to the url, we get Clue #2:

Clue 2

Again, very easy, and all that’s needed now is the first n digits of pi instead of e. The first 10 digit prime number found in consecutive digits of pi turns out to be 5926535897 (occurs very early!), leading to the final Clue #3.

Clue 3

I’ll leave you to do this one, and appreciate the punchline.

Code below the fold
Read the rest of this entry »

Written by parseval

March 6, 2009 at 11:40 pm

Posted in mathematics, programming

Tagged with

How the leopard may have really got her spots

with 3 comments

A mechanism of pattern formation in many animals, (including the zebra, leopard and the giraffe), was first suggested by Alan Turing in 1952.

To understand his mechanism, we’ll first look at a system where two reactive chemical species (morphogens) are present and they do not move about in space(ie, they don’t diffuse). If the concentration of the species are denoted by u and v, the rate at which the concentrations change will simply be the rate at which they react. That is,

\frac{\partial u}{\partial t} = f(u,v)

\frac{\partial v}{\partial t} = g(u,v)

where, f(u,v) and g(u,v) describe the rate law kinetics which the species obey. If we want to find the steady state concentrations of the species, we simply set the partial derivatives with time as 0, and solve for u and v.

What happens if f(u,v) and g(u,v) are non-linear equations with multiple solutions? The stability of each solution can be determined by a linear stability analysis, and the final stable solution will depend on the initial condition. Indeed, because of the homogeneity in the initial constraint, that the initial concentrations of the species are uniform everywhere and the they are spatially fixed, the final concentrations at steady state are again uniform in space.

Linearizing the above system of equations (ie, set u=u_{ss}+\tilde{u} and v=v_{ss}+\tilde{v}), one obtains

\left(\begin{array}{c} \frac{d \tilde{u}}{dt} \\ \frac{d \tilde{v}}{dt} \end{array}\right)=\left(\begin{array}{cc} f_u & f_v \\ g_u & g_v \end{array}\right) \left(\begin{array}{c} \tilde{u} \\ \tilde{v}\end{array}\right)

For the steady states to be stable, the condition is that the real parts of the eigenvalues of the jacobian are negative. For a two component system, this is true when the trace is negative and the determinant positive.

Now, what happens if we remove the spatial constraint and allow the chemical species to diffuse? The equations which describe the concentration of the species will be modified to account for diffusion. We’ll now have

\frac{\partial u}{\partial t} = D_u \nabla^2 u + f(u,v)

\frac{\partial v}{\partial t} = D_v \nabla^2 v + g(u,v)

where, D_u and D_v are the diffusivity of the chemical species. One of the observations which Turing made, was that under certain conditions, steady states which were initially spatially uniform could now show spatial variation due to diffusion. The natural question would be, under what conditions would we get these spatial variations? To answer that, we need to delve a bit into linear stability analysis.

However, before doing that, I’ll have to define what the boundary conditions and the geometry of the system is for the second set of equations. A very simple system would be a one-dimensional strip of space, where diffusion and reaction occurs. To make this system ‘self-contained’, we set the flux of the chemical species at the boundary to be 0.

Therefore, the system of equations are

\frac{\partial u}{\partial t} = D_u \frac{\partial^2 u}{\partial x^2} + f(u,v)

\frac{\partial v}{\partial t} = D_v \frac{\partial^2 v}{\partial x^2} + g(u,v)

Since we want to find out when the homogeneous solutions are unstable, we need to look at the linearized system in terms of the deviation variables from the steady state solutions, \tilde{u} = u - u_{ss} and \tilde{v} = v - v_{ss}. With these variables, the linearized equations are

\frac{\partial \tilde{u}}{\partial t} = D_u \frac{\partial^2 \tilde{u}}{\partial x^2} + f_u \tilde{u} + f_v\tilde{v}

\frac{\partial \tilde{v}}{\partial t} = D_v \frac{\partial^2 \tilde{v}}{\partial x^2} + g_u \tilde{u} + g_v\tilde{v}

The solutions to this equation can be shown to be of the form.

\tilde{u} = u^* e^{\sigma t} \sin{\alpha x}

\tilde{v} = v^* e^{\sigma t} \sin{\alpha x}

(Why? Substitute and check!)

Substituting to the equation, we get

u^* (\sigma +D_u\alpha^2 - f_u)  -f_v v^* = 0

v^* (\sigma +D_v\alpha^2 - g_v)  -g_u u^* = 0

Here’s the important argument. For a spatially varying pattern, we require non-zero solutions to u^* and v^*. The condition for that is

\left|\begin{array}{cc} \sigma +D_u\alpha^2 - f_u & -f_v \\ -g_u & \sigma +D_v\alpha^2 - g_v  \end{array}\right| = 0

This can be re-written as

det( \sigma I - A + \alpha^2 D) = 0


A= \left(\begin{array}{cc} f_u & f_v \\ g_u & g_v \end{array}\right)

D = \left(\begin{array}{cc} D_u & 0 \\ 0 & D_v \end{array}\right)

Now, the value of \sigma is the eigenvalue of the matrix R= (A - \alpha^2 D), and for the homogeneous solution to be unstable, one requires that the real part of the eigenvalues of R are positive for some \alpha^2

This condition can be simplified to,

\sigma^2 - Tr(R) \sigma + det(R) =0

For the eigenvalue to be positive, we require that det(R)<0 and Tr(R)>0 .

But, since Tr(R) = Tr(A - \alpha^2 D) = Tr(A) - \alpha^2 Tr(D) and Tr(A) is less than 0 because the homogeneous steady state is assumed to be stable, Tr(R) can never be positive. So, the criteria for instability of the homogeneous solution is that det(R)<0. One can also relate this to the ratio of the diffusivites of the two species u and v.

This means that under this condition, the homogeneous steady state solution is unstable, and the concentrations of u and v will increase till the non-linear terms we have neglected plays a role and causes spatial patterns.

I will edit and add pretty pictures displaying the simulation of this phenomenon once I figure out how to use the pdetool in MATLAB to solve non-linear parabolic partial differential equations in 2 variables.

The Chemical Basis of Morphogenisis, A.M Turing, Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, Vol. 237, No. 641. (Aug. 14, 1952),pp 37-72

Two-stage Turing model for generating pigment patterns on the leopard and the jaguar, Phys. Rev. E 74, 011914 (2006)

Mathematical Biology, JD Murray, 3rd edn, Springer (2002). chapter 2.4 &2.5

Written by parseval

April 14, 2008 at 4:28 pm

Nerd Sniping

with 5 comments

From xkcd,

Click to actually read it

So, how does one find the effective resistance?

The solution to this question, and many other network shapes, is present in Jozsef Cserti’s paper1.

Additionally, this website2 also clearly explains the technique to find the general equivalent resistance between any two nodes. The resistance between the two marked nodes in the xkcd question is,

\frac{1}{4 \pi^2} \int_{-\pi}^{\pi} \int_{-\pi}^{\pi} \left( \frac{1-\cos(2x+y)}{2-\cos{x} - \cos{y}} \right) dx dy

which turns out to be \frac{8-\pi}{2 \pi}. Check out Appendix A of Cserti’s paper for a technique to evaluate the above integral.

External Links
[1] – Application of the lattice Green’s function for calculating the resistance of an infinite network of resistors, Cserti József, American Journal of Physics, Volume 68, Issue 10, pp. 896-906 (2000). arXiv:cond-mat/9909120v4
[2] – Infinite 2D square grid of 1 ohm resistors

Written by parseval

January 3, 2008 at 12:02 am

Posted in humour, mathematics

Newton’s method and fractals

with one comment

Many of you would have come across the Newton-Raphson iterative technique of finding the root of a single non-linear equation, f(x)=0 in high-school.

The iterative method of solution is essentially given by,

x_{k+1}= x_{k} - \frac{f(x_k)}{f'(x_k)}, \quad k=0,1,2,....

Newton’s method can also be extended to find solutions to a set of simultaneous non-linear equations for functions of many variables.

For example, let’s say you need to find a solution to the equations,

\begin{cases} f_1(x_1,x_2,...,x_n)=0 \\ f_2(x_1,x_2,...,x_n)=0 \\ ... \\ f_n(x_1,x_2,...,x_n)=0 \end{cases}

with an initial guess, \bf{x} = (x^0_1,x^0_2,...,x^0_n)

The recursion which is used to find the solution is now given by

\bf{x}^{(k+1)} = \bf{x}^{(k)} - [\bf{J}(\bf{x}^{(k)}]^{-1} f(\bf{x}^{(k)}), \quad k=0,1,2,...

That is,

\left(\begin{array}{c} x^{(k+1)}_1 \\ x^{(k+1)}_2 \\ . \\x^{(k+1)}_n \end{array} \right) = \left( \begin{array}{c} x^{(k)}_1 \\ x^{(k)}_2 \\ . \\x^{(k)}_n \end{array} \right) - \left( \begin{array}{cccc} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & . & \frac{\partial f_1}{\partial x_n} \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & . & \frac{\partial f_2}{\partial x_n} \\ . & . & . & . \\ \frac{\partial f_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & . & \frac{\partial f_n}{\partial x_n} \end{array} \right)^{-1} \left(\begin{array}{c} f_1(\bf{x}^{(k)}) \\ f_2(\bf{x}^{(k)}) \\ . \\ f_n(\bf{x}^{(k)}) \end{array}\right)

Note that you’ll have to evaluate the Jacobian matrix at \bf{x} = \bf{x^{(k)}} during each iteration. For a proof of when the method converges, and some modifications to make it converge faster, check the references which are given at the end of this post.

There’s something interesting which happens when you apply this technique to an equation in the complex plane. For a complex number z=x+iy, any polynomial equation f(z) = 0 can be written as two equations in the variables x and y. For example, if f(z) = z^2 + 1 = 0, then,

(x+iy)^2 + 1 = 0

(x^2 -y^2 + 1) + i(2xy) = 0,

which means

f(x,y)=x^2-y^2+1=0 and g(x,y)=2xy=0

Also, note that for every (x,y) we have a corresponding point in the complex plane.

Now consider this question. Let’s say we have a polynomial equation z^n - 1 = 0. Given an arbitrary initial guess of the root to this equation

a) Will the algorithm converge to give a solution for that particular initial guess?
b) If so, how long does it take to converge to a solution for that particular initial guess?

An attempt to answer to these questions leads to fractal patterns. In the color scheme I’ve used, the bluer regions correspond to where the initial guesses converge faster. The yellow/red regions is where the initial guesses converge slowly or do not converge.







As you can see, while some initial guesses converge without much fuss (the basin regions), others can show remarkable behavior (guesses near the “boundary”) and can lead to a number of intricate patterns for various functions. I’ve linked to some websites which contain a gallery of such images for various functions. Some of them are truly remarkable!

References and External Links
A Class of Methods for Solving Nonlinear Simultaneous Equations, C. G. Broyden, Mathematics of Computation, Vol. 19, No. 92 (Oct., 1965), pp. 577-593
– An Introduction to Numerical Analysis, Endre Suli and David Mayers, Cambridge University Press,Pg 104-125
Fractals from Newton’s Method
Newton Basins
Newton’s method and Newton basin fractals
Newton’s Method gallery
Another Newton’s Method Gallery

Written by parseval

September 20, 2007 at 1:23 pm

Posted in mathematics

Ain’t so convoluted

with 3 comments

What’s a convolution, you ask? If you have two functions f(x) and g(x), the convolution of f and g, is defined as

f \otimes g = \int\limits_{-\infty}^{\infty} f(u) g(x-u) du

where u is the dummy variable of integration.

The notation f \ast g is also used to denote the convolution of f and g, although I’ll be using the former notation. Note that the convolution of f and g is itself a function of x. That is,

h(x) = f \otimes g

Here’s one way to visualize the convolution. If you have two functions f and g with respect to the dummy variable, what the convolution actually does is,

STEP 1: “Flips” one of the functions about the y-axis. For example, g(u) to g(-u)

STEP 2: Shifts it by an amount x. ie, to g(x-u)

STEP 3: “Slides” g(x-u) along the u-axis by keeping u fixed, and allowing x to vary all the way from -\infty to \infty. The value of the convolution at some point x_0, is h(x_0) = \int\limits_{-\infty}^{\infty} f(u)g(x_0 - u) du

So, if you consider the convolution integral, the value of the product h(u) g(x_0-u), and hence h(x_0) is zero when the two functions do not intersect. However, when the two functions do intersect, the value of the convolution at that point will be the integral of the product over the entire overlapping region (where the product is non-zero), and this value is simply the area of the overlapping region.

You can see detailed animated illustrations of this idea at Wolfram’s MathWorld.

Another useful way of thinking about the convolution of two functions, is by the concept of a functional. We can consider h as a functional of the function f. That is, for every given function f, there will be a corresponding value of h. So, to calculate the value of h(x) at some point x_0, we still need to know the entire function f.

There are useful properties associated with convolution. For example, convolutions are
(i) Commutative

f \otimes g = g \otimes f

So, it doesn’t matter which function you flip.

(ii) Associative

f \otimes (g \otimes h) = (f \otimes g) \otimes h

(iii) Distributive

f \otimes (g+h) = f \otimes g + f \otimes h

Convolution theorem:
The convolution theorem is immensely useful to calculate fourier transforms and find fourier pairs. It states that, if F_1(x) \rightleftharpoons \phi_1(p) and F_2(x) \rightleftharpoons \phi_2(p),

F_1(x) \otimes F_2(x) \rightleftharpoons \phi_1(p) \cdot \phi_2(p)

It’s an interesting result that the fourier transform of the convolution of two functions, is the product of the corresponding fourier transforms of the individual functions. This means that a convolution in the normal domain, becomes a product in the fourier domain and vice-versa

It’s actually easy to prove this useful result. For, if the convolution of two functions is given as
h(x) = \int\limits_{-\infty}^{\infty}F_1(u)F_2(x-u)du

Then the fourier transform of h(x) is

\int\limits_{-\infty}^{\infty} h(x) e^{2 \pi i p x} dx

= \int\limits_{-\infty}^{\infty} \left( \int\limits_{-\infty}^{\infty} F_1(u)F_2(x-u) du\right) e^{2 \pi i x p} dx

= \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} F_1(u) F_2(x-u) e^{2 \pi i x p} du dx

Since x is a dummy variable, set x-u=z, while holding u constant, so that the integral reduces to

\int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} F_1(u) F_2(z) e^{2 \pi i p (z+u)} du dz

= \left(\int\limits_{-\infty}^{\infty} F_1(u) e^{2 \pi i p u} du\right) \cdot \left( \int\limits_{-\infty}^{\infty} F_2(z) e^{2 \pi i z p} dz \right)

= \phi_1(p) \cdot \phi_2(p)

We can apply this theorem to find various convolutions. Here’s an example.

(i) Convolution with a delta function

Let f(x) \rightleftharpoons \phi(p)

We’re now interested in the convolution of f with a shifted delta function, f(x) \otimes \delta(x-a)
To find this, let’s apply the convolution theorem. We know that f(x) \rightleftharpoons \phi(p) and \delta(x-a) \rightleftharpoons e^{- 2 \pi i p a}. Therefore f(x) \otimes \delta(x-a) \rightleftharpoons \phi(p) \cdot e^{- 2 \pi i p a}

And taking the inverse fourier transform of \phi(p) \cdot e^{- 2 \pi i p a}, we find that f(x) \otimes \delta(x-a) is f(x-a). Notice that f(x) is shifted by an amount a.

Finally, I’ll discuss one more property I’ll need to use before I move on to a part of what I’m currently working on, which is the application of fourier analysis in interferometry and diffraction. This property is the derivative of fourier transforms.

If, F(x) \rightleftharpoons \phi(p),
\frac{dF}{dx} \rightleftharpoons -2 \pi i p \phi(p)

and therefore,

\frac{d^n F}{dx^n} \rightleftharpoons \left(-2 \pi i p \right)^n \phi(p)

which you can easily check this, by use of the Leibniz rule. This property of fourier transforms is useful in solving linear PDE’s.

Written by parseval

June 30, 2007 at 8:18 pm

Posted in mathematics

FT of some common functions

with one comment

In this post, let’s look at the fourier transform of some functions which are quite useful.

(i) The rectangle function
We’ll start with the rectangular function, also called the box function. It’s defined as

f_a(x) = \begin{cases} 0, \quad |x| > a/2 \\ 1, \quad |x| \leq a/2 \end{cases}

The Box hat function

Now, consider the fourier pair of f(x). We have,

\phi(p) = \int\limits_{-\infty}^{\infty} f(x) e^{2 \pi i p x}dx

= \int\limits_{-a/2}^{a/2} e^{2 \pi i p x}dx

= \left(e^{2 \pi i p a/2} - e^{-2 \pi i p a/2}\right)/ \left(2 \pi i p \right)

= a \sin(\pi p a)/(\pi p a)

= a \text{sinc} (\pi p a)

So, the fourier transform of the boxcar function is the sinc function!

Sinc function

I’ll revist this fourier pair again, while discussing the wave theory of light. In fact, we can use this fourier pair to show that the interference pattern we get in a double slit experiment is infact the sinc function!

(ii) The Gaussian function

Next, we’ll look at the gaussian function. The gaussian function has the interesting property that it’s fourier pair is also a gaussian function! Consider,

f(x) = e^{-x^2/a^2}


Let’s caculate the fourier pair, \phi(p).

\phi(p) = \int\limits_{-\infty}^{\infty} e^{\frac{-x^2}{a^2}} e^{2 \pi i p x} dx
= \phi(p) = \int\limits_{-\infty}^{\infty} e^{-\left(x/a - a \pi i p\right)^2} e^{-pi^2 p^2 a^2} dx

To evaluate this integral, use the substitution,
x/a - a \pi i p=r

After evaluating the integral and substituting the limits, the expression for \phi(p) is obtained as,

\phi(p) = a \sqrt{\pi} e^{-\pi^2 p^2 a^2}


which is also a gaussian.

Also, if you try plotting the fourier pairs for different values of a (and hence, different widths of the gaussian), you’ll notice that the wider the gaussian in x-space, the narrower it is in p-space (ie, the transform space), and vice versa.

(iii) The delta function

The dirac delta function (although, not strictly a function), can be represented as
\delta(x) = \begin{cases} 0, \quad x \neq 0 \\ \infty, \quad x = 0 \end{cases}

Now let’s apply the fourier transform to the delta function. We get,

\phi(p) = \int\limits_{-\infty}^{\infty} \delta(x) e^{2 \pi i p x} dx

and by the property of the delta function, this is,

= e^{2 \pi i p 0}

= 1

Therefore, we find that the fourier pair of the delta function is unity. That is \delta(x) \rightleftharpoons 1

Also, notice that
\delta(x-a) \rightleftharpoons e^{2 \pi p a}
\delta(x+a) \rightleftharpoons e^{- 2 \pi p a}

and hence,
\delta(x+a) + \delta(x-a) \rightleftharpoons 2 \cos(2 \pi p a)

(iv) The Shah function

The shah function, also known as a Dirac comb, is an infinite combination of evenly spaced dirac functions.

f_a(x)= \sum_{n=-\infty}^{\infty} \delta(x-an)

The fourier transform of the shah function is also another shah function, with a period of 1/a. You’ll find that the shah function is quite invaluable in convolutions, where it’s role is to create infinte “copies” of the original function, with period equal to the spacing between the teeth of the comb.

In my next post, I’ll explain more about convolution, especially the convolution theorem which is a real time saver in performing transforms, and other theorems relating to fourier transforms

Note: If you find any errors, please do inform me, and I’ll correct them. Also, click on the thumbnail images to get a detailed graph

Written by parseval

June 22, 2007 at 2:26 pm

Posted in mathematics

Fourier transforms for the practical person

with 3 comments

In this series of posts, I plan to outline some basic ideas which I’ve learnt on the theory of Fourier transforms, and it’s practical applications in a non-rigorous manner. Once I’ve laid out the basics, I’ll then show you some interesting stuff from what I’m currently working on.

Let’s start with Fourier series. It’s actually a remarkable fact, that we can express any arbitrary periodic function, simply as the sum of the ordinary sine and cosine functions we’ve all studied at high school. If F(t) is a periodic function, then we have

F(t) = \sum_{n=-\infty}^{\infty} A_n \cos(2 \pi n \nu_0 t) + B_n \sin(2 \pi n \nu_0 t)

The Fourier transform is an extension of this, as the period of the function approaches infinity, and the gap between successive harmonics approaches 0. So, in some sense, the fourier transform decomposes a function into it’s frequency components.

For a non-periodic function F(t) which satisfies certain conditions, there are many conventions of describing the fourier transform. Following one such convention which is widely used, the forward fourier transform is

F(t)=\int\limits_{-\infty}^{\infty} \phi(\nu) e^{-2 \pi i \nu t} d\nu

While, the inverse fourier transform is

\phi(\nu)=\int\limits_{-\infty}^{\infty} F(t) e^{+2 \pi i \nu t} dt

Notice that, if F(t) is a continuous time signal, then it’s transformed into the frequency domain by the forward transform. One of the properties of the fourier transforms is that, \phi(\nu) and F(t) are transforms of each other, and form a fourier pair, and are represented by F(t) \rightleftharpoons \phi(\nu).

This means that, if f(x) \rightleftharpoons \phi(p) , then

f(x)=\int\limits_{-\infty}^{\infty} \left( \int\limits_{-\infty}^{\infty} f(x) e^{2 \pi i x p}dx \right) e^{-2 \pi i x p} dp

if f(x) isn’t discontinuous. If it is discontinuous, then the value at that point will be the average of the value around the discontinuity. So, we can simplify our terminology and say that the fourier transform of \phi(p) is f(x) and vice versa.

In the next post, I’ll look at the fourier transform of some useful functions, but before that, there’s one more nice result. For an Electromagnetic wave, or a signal in a wire, the fourier transform of the voltage can be complex. However the conjugate product \phi(\nu) \phi^{\ast}(\nu)=|\phi(\nu)|^2, is real and is proportional to the power density (or, power per unit frequency). This is know as the spectral power density.

Written by parseval

June 16, 2007 at 8:35 am

Posted in mathematics