Pi Day Post

March 14th marks the annual celebration of the beloved circumference to diameter ratio, π! Unless you're British, in which case, see you on July 22.




The essence of today's celebration is in the American style date, written 3-14, which lists the first three numbers in the mathematical constant. Last year, in 2015, the pi day was especially accurate, since the date was 3-14-15, the first five digits of the constant. Even better were the folks who celebrated the never-to-be-seen-again pi second, at precisely 3-14-15 9:26:53. But I object to that exclusivity, because today's pi day, in 2016, is just as special!

MATLAB by default actually uses an approximation for pi, and as you probably guessed by now, it's: $$3.1416$$
This corresponds to today's date, without the need to get specific about seconds! It will be pi day all day!

So, let's make some scripts in MATLAB that approximate PI. This first one is a bit complicated, and isn't that great, but I love its implementation of prime numbers. This sequence says (I haven't found an official formula for it; I am too bad at math to figure one out for myself!) for an infinite series of prime numbers (it is actually proven that an infinite amount exist! Thanks, Euclid!):

• reciprocalize each term (1 / prime)
• if the prime in the term is one more than a factor of 4, (prime % 4 == 1), add it to 1
• otherwise, subtract it from 1
• multiply every term

This should yield a product that is approximately equivalent to $$\prod = \frac{2}{\pi}$$

However, the approximation isn't that great, and gets closest at a maximum prime-value of 475.

	max = input('Maximum size of primes: ');
	h = zeros(1, max);
	%load array with harmonic series
	for i = 1 : max
	h(i) = i^(-1);
	end
	sum = 0;
	%add each term squared
	for i = 1 : max
	sum = sum + h(i)^2;
	end
	%for a finite number of terms, the sum should be approximately pi^2 / 6
	sum = sqrt(6 * sum);
	fprintf('With %.f terms in the sum, pi estimates to be %f!\n', max, sum);

view raw pi_day1.m hosted with ❤ by GitHub

The second algorithm is actually much more effective and simple! We take the standard Harmonic Series and:

• square each term
• sum each term 

This sum gives us $$\sum_{n=1}^{+\infty} \frac{1}{n^2} = \frac{\pi^2}{6}$$

	max = input('Maximum size of primes: ');
	h = zeros(1, max);
	%load array with harmonic series
	for i = 1 : max
	h(i) = i^(-1);
	end
	sum = 0;
	%add each term squared
	for i = 1 : max
	sum = sum + h(i)^2;
	end
	%for a finite number of terms, the sum should be approximately pi^2 / 6
	sum = sqrt(6 * sum);
	fprintf('With %.f terms in the sum, pi estimates to be %f!\n', max, sum);

view raw pi_day2.m hosted with ❤ by GitHub