Create plots of the convergence of Madhava's series and Gregory's series

plot_series_convergence.py

#!/usr/bin/env python3
"""
plot_series_convergence.py

Plot Gregory's series and Madhava's series
to demonstrate their convergence speed.

Here we refer to Madhava's series as the rapidly converging series
    pi = sqrt(12) * (1 - 1/(3 * 3)  + 1/(3 * 5^2) - 1/(3 * 7^3))
which is described here:
    https://en.wikipedia.org/wiki/Madhava_of_Sangamagrama#The_value_of_%CF%80_(pi)
Gregory's series was also discovered by Madhava (and Leibniz): it converges
much more slowly:
    https://en.wikipedia.org/wiki/Gregory%27s_series
"""
import matplotlib.pyplot as plt
import math


TEST_GREGORYS = False
PLOT_GREGORYS = True
TEST_MADHAVAS = False
PLOT_MADHAVAS = True
PLOT_COMPARISON = True


def madhavas_series(n):
    """
    Calculate Madhava's series to n terms:
        pi = sqrt(12) * (1 - 1/(3 * 3)  + 1/(3 * 5^2) - 1/(3 * 7^3))
    This is the rapidly-converging series described here:
        https://en.wikipedia.org/wiki/Madhava_of_Sangamagrama#The_value_of_%CF%80_(pi)
    """
    result = 0
    sign = 1
    next_odd = 1
    exponent = 0
    for i in range(n):
        result = result + sign * (1/(3**exponent * next_odd))
        exponent += 1
        next_odd += 2
        sign = -sign
    return math.sqrt(12) * result


def gregorys_series(n):
    """
    Calculate Gregory's series to n terms:
        pi/4  =  1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 ...
    https://en.wikipedia.org/wiki/Gregory%27s_series
    """
    result = 0
    next_odd = 1
    sign = 1
    for i in range(n):
        result += sign / next_odd
        next_odd += 2
        sign = -sign
    return 4 * result


def test_gregorys_series():
    """
    Demonstrate the convergence of Gregory's series
    non-visually.
    """
    print('-------------------------------------------')
    for l in range(1, 6):
        n = 10**l  # Increment by powers of ten
        result = gregorys_series(n)
        error = result - math.pi
        print('Gregory\'s series: {} terms evaluated'.format(n))
        print('Result: %.18f' % result)
        print('Error:  %.18f' % error)
        print('-------------------------------------------')


def test_madhavas_series():
    """
    Demonstrate the convergence of Madhava's series
    non-visually.
    """
    print('-------------------------------------------')
    for l in range(1, 6):
        n = 10**l  # Increment by powers of ten
        result = madhavas_series(n)
        error = result - math.pi
        print('Madhava\'s series: {} terms evaluated'.format(n))
        print('Result: %.18f' % result)
        print('Error:  %.18f' % error)
        print('-------------------------------------------')


def pi_plot(n, series, series_name):
    """
    Plots the first n iterations of the series specified to show
    the rate a which it converges to π
    Also, its a pun.
    """
    fig = plt.figure()
    ax = fig.add_subplot()
    iter_labels = [i for i in range(n)]
    results = []
    for i in range(n):
        results.append(series(i))
    ax.bar(iter_labels, results, align='center', color='b', alpha=0.85)
    ax.axhline(y=math.pi, color='r', linestyle='-', label='pi')
    ax.set_xlabel('Iterations')
    ax.set_title(series_name)
    ax.legend()
    plt.show()


def pi_plot_compare(n, s1, s2, s1name, s2name):
    """
    Compare the convergence plots of two series on the same figure.
    """
    width = 0.25
    fig = plt.figure()
    ax = fig.add_subplot()
    s1_iter = [i for i in range(n)]
    s1_results = []
    for i in range(n):
        s1_results.append(s1(i))
    ax.bar(s1_iter, s1_results, width=width, color='g', alpha=0.85, label=s1name)
    s2_iter = [i+width for i in range(n)]
    s2_results = []
    for i in range(n):
        s2_results.append(s2(i))
    ax.bar(s2_iter, s2_results, width=width, color='b', alpha=0.85, label=s2name)
    ax.axhline(y=math.pi, color='r', linestyle='-', label='pi')
    ax.set_xlabel('Iterations')
    ax.set_title(s1name + ' vs. ' + s2name)
    ax.legend()
    plt.show()


def main():
    """
    Run each function specified in the global flags.
    """
    if TEST_MADHAVAS:
        test_madhavas_series()
    if TEST_GREGORYS:
        test_gregorys_series()
    if PLOT_MADHAVAS:
        pi_plot(30, madhavas_series, 'Madhava\'s series')
    if PLOT_GREGORYS:
        pi_plot(40, gregorys_series, 'Gregory\'s series')
    if PLOT_COMPARISON:
        pi_plot_compare(30, madhavas_series, gregorys_series,
                        'Madhava\'s series',
                        'Gregory\'s series')


if __name__ == "__main__":
    main()