Last active
March 1, 2018 16:02
-
-
Save person142/9789b09c599bda1c2377d496787d4ea0 to your computer and use it in GitHub Desktop.
Make sure Halley's method gets cubic convergence
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
import mpmath | |
from scipy.optimize import newton | |
import matplotlib.pyplot as plt | |
points = [] | |
def f(x): | |
points.append(x) | |
return mpmath.sin(x) | |
def fprime(x): | |
return mpmath.cos(x) | |
def fprime2(x): | |
return -mpmath.sin(x) | |
def main(): | |
with mpmath.workdps(4000): | |
tol = mpmath.mpf('1e-4000') | |
try: | |
root = newton(f, mpmath.mpf(3.5), fprime=fprime, fprime2=fprime2, | |
tol=tol, maxiter=100) | |
except RuntimeError: | |
pass | |
digits = [] # roughly the number of correct digits | |
for p in points: | |
digits.append(abs(mpmath.log10(abs(p - mpmath.pi)))) | |
ratios = [] # they should all be close to 3 | |
for i in range(len(digits) - 1): | |
ratios.append(digits[i+1]/digits[i]) | |
plt.plot(ratios) | |
plt.show() | |
if __name__ == '__main__': | |
main() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment