Created
November 10, 2019 09:32
-
-
Save eindiran/e110720fc16fe07ca11dae5b8c371308 to your computer and use it in GitHub Desktop.
Reimplementation of plot_series_convergence.py using mpmath library for high precision floating-point arithmetic.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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 | |
from mpmath import * | |
mp.dps = 75 # 75 decimals of precision | |
TEST_GREGORYS = True | |
PLOT_GREGORYS = True | |
TEST_MADHAVAS = True | |
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 = fadd(result, fmul(sign, fdiv(1, fmul(3**exponent, next_odd)))) | |
exponent = fadd(exponent, 1) | |
next_odd = fadd(next_odd, 2) | |
sign = fneg(sign) | |
return fmul(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 = fadd(result, fmul(sign, fdiv(1, next_odd))) | |
next_odd = fadd(next_odd, 2) | |
sign = fneg(sign) | |
return fmul(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 = fsub(result, pi) | |
print('Gregory\'s series: {} terms evaluated'.format(n)) | |
print('Result: %s' % nstr(result, n=25)) | |
print('Error: %s' % nstr(error, n=25)) | |
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 = fsub(result, pi) | |
print('Madhava\'s series: {} terms evaluated'.format(n)) | |
print('Result: %s' % nstr(result, n=25)) | |
print('Error: %s' % nstr(error, n=25)) | |
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=pi, color='r', linestyle='-', label='pi') | |
ax.set_xlabel('Iterations') | |
ax.set_title(series_name) | |
ax.set_ylim(2.6, ax.get_ylim()[1]) | |
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=pi, color='r', linestyle='-', label='pi') | |
ax.set_xlabel('Iterations') | |
ax.set_title(s1name + ' vs. ' + s2name) | |
ax.set_ylim(2.6, ax.get_ylim()[1]) | |
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(100, 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() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment