Last active
November 12, 2021 15:14
-
-
Save pierre-haessig/42972f91dae1edd34f2df86462f6280a to your computer and use it in GitHub Desktop.
Numerical implementation of the formula for the partial sum of a geometric series when ratio equals one. See https://math.stackexchange.com/questions/4288069/formula-for-the-partial-sum-of-a-geometric-series-when-ratio-equals-one
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
# Comparison of three numerical implementation of the formula for the partial sum of a geometric series | |
# when ratio equals (or is close to) one. | |
# See https://math.stackexchange.com/questions/4288069/formula-for-the-partial-sum-of-a-geometric-series-when-ratio-equals-one | |
from math import expm1, log1p | |
def SN_def(a,N): | |
return sum(a**n for n in range(N)) | |
def SN_form(a,N): | |
return (1-a**N)/(1-a) | |
def SN_form2(eps,N): | |
return -expm1(N*log1p(-eps))/eps | |
def SN_approx(eps,N): | |
return N*(1-(N-1)/2*eps) | |
# Compare the implementations: | |
N = 100 | |
eps_list = [10**(-i) for i in [1, 3, 5, 7, 9, 11, 13, 15, 17, 19]] | |
SN_fmt = lambda x: '{:20.15f}'.format(x) | |
txt_fmt = lambda s: '{:>20s}'.format(s) | |
sep = ', ' | |
print(' '*7, end=sep) | |
print(txt_fmt('Definition with sum'), end=sep) | |
print(txt_fmt('Formula with expm1'), end=sep) | |
print(txt_fmt('Approx for small ε'), end=sep) | |
print(txt_fmt('Classical formula')) | |
for eps in eps_list: | |
a = 1-eps | |
print('ε={:.0e}'.format(eps), end=sep) | |
print(SN_fmt(SN_def(a,N)), end=sep) | |
print(SN_fmt(SN_form2(eps,N)), end=sep) | |
print(SN_fmt(SN_approx(eps,N)), end=sep) | |
try: | |
print(SN_fmt(SN_form(a,N))) | |
except ZeroDivisionError: | |
print(txt_fmt('zero div. error')) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment