A simple convergence comparison between Leibniz's and mine pi formula

Author: Kazi Abu Rousan

Here there is $\pi$, there is circle.

Today's blog will be a bit different. This will not discuss any formula or any proof, rather it will just contain a python program to compare a $\pi$ formula given by Leibniz and a one discovered by me (blush).

How does Leibniz formula looks like?

The Leibniz formula for $\pi$ is named after Gottfried Leibniz( published in 1676). Sometimes it is also called Leibniz-Madhava series as this is actually a special case of a more general formula discovered by Madhava of Sangamagrama in 14th century.

The formula is:

$\frac{\pi}{4}$ = $1 - \frac{1}{3}$ + $\frac{1}{5}$ - $\frac{1}{7}$ + $\frac{1}{9}$ -$\cdots$ + $\frac{1}{4n+1}$ - $\frac{1}{4n+3}$ + $\cdots$

If you want a visual proof, read my article : Article- click

or this video: Geometric proof

So, how can you find $\pi$ using this?, Just calculate each fractions and add.

As an example, $4, 4(1-\frac{1}{3}) = 2.6667$, $4(1-\frac{1}{3}$ + $\frac{1}{5}) = 3.4664$, $4(1-\frac{1}{3}$ + $\frac{1}{5}$ - $\frac{1}{7}) =2.895238 $, just like this.

But it is nowhere close to the value of $\pi = 3.141592653\cdots$

The conclusion is that, it converges very slowly. Let's ask a computer to find the value of $\pi$ using the series.

pi_val = 3.1415926535897932#right value upto 17 decimals
calculated_pi_from_series = 0
n = 1000
for i in range(0,n):
    calculated_pi_from_series = calculated_pi_from_series + 1/(4*i+1) - 1/(4*i+3)
calculated_pi_from_series = 4*calculated_pi_from_series
print(calculated_pi_from_series)
print("error = ",pi_val - calculated_pi_from_series)
#output: 3.1410926536210413
#output: error =  0.0004999999687518297

For a output of $3.141592153589724$ and error $= 5.000000689037165e-07$, we need $n$ = $10,00,000$. As you can see slow this series is. But even then it has great mathematical significance. It even has a beautiful geometric relationship with Fermat's two square theorem and prime numbers.

How does my formula looks like?

As it is my own discovery, it doesn't have any name but you may call it Rousan's Formula (🤣🤣). This one is almost identical to Leibniz's formula. It is so as it also uses 2d complex numbers as it's basis just like Leibniz's formula.

The main difference is that it mine uses Eisenstein Integers, i.e., triangular grid while Leibniz's uses Gaussian Integers, i.e., regular square grid.

This formula can be written as :

$\frac{\pi}{3\sqrt{3}}$ = $1 - \frac{1}{2}$ + $\frac{1}{4}$ - $\frac{1}{5}$ + $\frac{1}{7}$ - $\cdots$ + $\frac{1}{3n+1}$ - $\frac{1}{3n+2}$ + $\cdots$

Almost same right?, now let's use python to see how fast it converges.

import math as mp
pi_val = 3.1415926535897932#right value upto 17 decimals
calculated_pi_from_series = 0
n = 1000
for i in range(0,n):
    calculated_pi_from_series = calculated_pi_from_series + 1/(3*i+1) - 1/(3*i+2)
calculated_pi_from_series = 3*calculated_pi_from_series*mp.sqrt(3)
print(calculated_pi_from_series)
print("error = ",pi_val - calculated_pi_from_series)
#output: 3.141015303363362
#output: error =  0.0005773502264312391

For $n = 10,00,000$, we get $3.1415920762392955$ with an error of $5.773504976325228e-07$. So, This shows mine converges a bit more slowly.

A direct comparison and reason

So, why mine converge more slowly?, well it is hidden in it's shape. Mine uses triangular grid and leibniz's uses squares. This is the main reason. But fear not, I have actually found a more generalized version of fermat's two square theorem, by using it i can get another similar series which converges faster than leibniz. Maybe I will discuss that in some later blog.

For now, let's try to plot both series's error. To we will plot the number of terms taken in x axis and error in y axis.

import matplotlib.pyplot as plt
import math as mp
pi_val = 3.1415926535897932
sqrt_3_3 = 3*mp.sqrt(3)

def leib_err(n):
    result = 0
    for i in range(n):
        result = result + 1/(4*i+1) - 1/(4*i+3)
    result = 4*result
    error = pi_val - result
    return result

def rou_err(n):
    result = 0
    for i in range(n):
        result = result + 1/(3*i+1) - 1/(3*i+2)
    result = result*sqrt_3_3
    error = pi_val - result
    return result
n_vals = [i for i in range(10_000+1)]
leibniz = [leib_err(i) for i in n_vals]
# print(leibniz)
rousan  =  [rou_err(i) for i in n_vals]
#
plt.plot(n_vals,leibniz,c='red',label='Leibniz Error')
plt.plot(n_vals,rousan,c='black',label='My formula Error')
plt.legend(loc='best',prop={'size':7})
plt.ylim(3.1375,3.1425)
plt.title("Error vs No. of terms for Leibniz and my $\\pi$ formula.")
plt.show()

The output is:

Error vs no. of terms for Leibniz and my pi formula

Error vs no. of terms for Leibniz and my pi formula - 2

This is all for today, see you guys next time with something new.

Your Personal Mandelbrot Set

Author: Kazi Abu Rousan

Bottomless wonders spring from simple rules, which are repeated without end.
Benoit Mandelbrot

Today, we will be discussing the idea for making a simple Mandelbrot Set using python's Matplotlib. This blog will not show you some crazy color scheme or such. But rather the most simple thing you can make from a simple code. So, Let's begin our discussion by understanding what is a Mandelbrot Set.

What is Mandelbrot Set?

In the most simple terms, The Mandelbrot set is the set of all complex numbers c for which the function $f_c(z) = z^2 + c$ doesn't goes to infinity when iterated from $z=0$.

Let's try to understand it using a example.

Suppose, we have a number $c = 1$, then to $ f_c(z)=f_1(z)= z^2 + 1 $.

Then, we first take $z_1 = f_1(z=0)=1$, after this we again calculate $z_2 = f_1(f_1(0)) = 1^2 + 1 = 2$, and so on. This goes on, i.e, we are getting new z values from the old one using $z_{n+1} = {z_n }^2 + 1$. This generates the series $0,1,2,5,26,\cdots$. As you can see, each element of this series goes to infinity, so the number $1$ is not an element of mandelbrot set.

Similarly, if we take $c=-1$, then the series is $0,-1,0.-1,0,\cdots$. This is a oscillating series and hence it will never grow towards infinity. Hence, the number $-1$ is an element of this set.

I hope this idea is now clear. It should be clear that $c$ can be any imaginary number. There is no restriction that says that $c$ must be real. So, is the number $c=i$ inside the set? Try it!!

Check if a Number is inside Mandelbrot set or not

Let's take a point $c = x+iy = P(x,y)$. So, the first point is $z_0 = 0+i\cdot 0 \to x = 0$ and $y=0$, i.e., we start the iteration from $f_c{z=0}$. Now, if at some point of iteration $|z_n|\geq 2$, then the corresponding $c$ value will not be inside the set. But if $|z|<2$, then corresponding c is inside the set.

Try this idea by hand for new numbers.

Now, let's check few numbers using a simple python code.

import numpy as np
def in_mandel(c,max_ita=1000):
    z = 0
    for i in range(max_ita):
        z = z**2 + c
        if abs(z)>=2:
            print(c," is not inside mandelbrot set upto a check of ",(i+1))
            break
        if i==max_ita-1 and abs(z)<2:
            print(c," is inside mandelbrot set upto a check of ",max_ita)
c1 = 1
in_mandel(c1)
#Output = 1  is not inside mandelbrot set
c2 = -1
in_mandel(c2)
#Output = -1  is inside mandelbrot set upto a check of  1000
c3 = 0.5+0.4j
in_mandel(c3)
#output = (0.5+0.4j)  is not inside mandelbrot set

You see how easy the code is!!, you can try out this geogebra version of you want - Check it here: Geogebra

Plot Mandelbrot set

Now, Let's plot one. We will be using a simple rule.

If any number c is inside the set, we will mark that coordinate in the graph.
If it is not there, we will just leave it.

import matplotlib.pyplot as plt
import numpy as np
max_ita = 1000 # number of interations
div = 1000 # divisions
c = np.linspace(-2,2,div)[:,np.newaxis] + 1j*np.linspace(-2,2,div)[np.newaxis,:]
z = 0 
for n in range(0,max_ita):
      z = z**2 + c
p = c[abs(z) < 2]#omit z which diverge 
plt.scatter(p.real, p.imag, color = "black" )
plt.xlabel("Real")
plt.grid(color='purple',linestyle='--')
plt.ylabel("i (imaginary)")
plt.xlim(-2,2)
plt.ylim(-1.5,1.5)
plt.show()

The output is:

This is not as good as the plots you normally see in internet like the images given below. But this one serves the purpose to show the shape of the set.

Let's try to make a bit fancier one.

The code is given here:

import matplotlib.pyplot as plt
import numpy as np
max_ita = 1000
div = 1000
c = np.linspace(-2,2,div)[:,np.newaxis] + 1j*np.linspace(-2,2,div)[np.newaxis,:] 
ones = np.ones(np.shape(c), np.int)
color = ones * max_ita + 5
z = 0
for n in range(0,max_ita):
      z = z**2 + c
      diverged = np.abs(z)>2
      color[diverged] = np.minimum(color[diverged], ones[diverged]*n)
plt.contourf(c.real, c.imag, color)
plt.xlabel("Real($c$)")
plt.ylabel("Imag($c$)")
plt.xlim(-2,2)
plt.ylim(-1.5,1.5)
plt.show()

Here is the output:

Beautiful isn't it?, This set has many mystery. Like you can find the value of $\pi$ and Fibonacci number from it. It is also hidden in Bifurcation diagram and hence related to chaos theory.

But we will live all these mysteries for some later time. This is all for today. I hope you have enjoyed this and if you have, then don't forget to share.