Make sure Halley's method gets cubic convergence

test_halley_convergence.py

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()