This file has many functions that can iterate and plots the final result for convergence, precision is also a input

CompareNumericalConvergence.py

import numpy as np
import sympy as sp
from matplotlib import pyplot as plt
import mpmath as mp
mp.mp.dps = 4000
print(mp.mp)
sqrt = mp.sqrt


def xsqrta(x, a):
    # return x * sqrt(a)
    return (1+a*x**2)/2

def bab(x, a):
    return (x+1/(x*a))/2

def bab2(x, a):
    x = bab(x, a)
    return bab(x, a)

def stac(x, a):
    return x*(3-a*x**2)/2

def g11(x, a):
    """f = x - 1/(a*x)"""
    return 2*x/(1+a*x**2)

def g120(x, a):
    """f = x - 1/(a*x)"""
    inner = (3+a*x**2)**2 - 12
    if inner < 0:
        msg = "Problem in g120:\n"
        msg += "    inner = %.4f\n" % inner
        msg += "        x = %.4f\n" % x
        msg += "        a = %.4f\n" % a
        msg += "   a*x**2 = %.4f\n" % (a*x**2)
        raise ValueError(msg)
    raiz = sqrt(inner)
    g = (x/2)*(3 + a*x**2 - raiz)
    return g

def g121(x, a):
    """f = x - 1/(a*x)"""
    # raiz = sqrt( (3+a*x**2)**2 - 12)
    raiz = 4*xsqrta(x, a)-2
    g = (x/2)*(3 + a*x**2 - raiz)
    return g

def g122(x, a):
    """f = x - 1/(a*x)"""
    # raiz = sqrt( (3+a*x**2)**2 - 12)
    raiz = 6*xsqrta(x, a) - (3+a*x**2)
    g = (x/2)*(3 + a*x**2 - raiz)
    return g

def g123(x, a):
    """f = x - 1/(a*x)"""
    # raiz = sqrt( (3+a*x**2)**2 - 12)
    raiz = 3*xsqrta(x, a)*(5+a*x**2) - 10*a*x**2 - 6
    g = (x/2)*(3 + a*x**2 - raiz)
    return g

def g21(x, a):
    """f = x^4 - 1/a^2"""
    return (3*x/4) + 1/(4*a**2*x**3)

def g220(x, a):
    """f = x^4 - 1/a^2"""
    inner = 2*(3 - a**2 * x**4)
    if inner < 0:
        msg = "Problem in g220:\n"
        msg += "    inner = %.4f\n" % inner
        msg += "        x = %.4f\n" % x
        msg += "        a = %.4f\n" % a
        msg += "     x**4 = %.4f\n" % (x**4)
        msg += "     a**2 = %.4f\n" % (a**2)
        raise ValueError(msg)
    raiz = sqrt(inner)
    return (2*x/3) + raiz/(6*a*x)

def g221(x, a):
    """f = x^4 - 1/a^2"""
    # raiz = sqrt(2*(3 - a**2 * x**4))
    raiz = 4-2*xsqrta(x, a)
    return (2*x/3) + raiz/(6*a*x)

def g222(x, a):
    """f = x^4 - 1/a^2"""
    # raiz = sqrt(2*(3 - a**2 * x**4))
    raiz = 6*xsqrta(x, a) - 4*a*x**2
    return (2*x/3) + raiz/(6*a*x)

def g223(x, a):
    """f = x^4 - 1/a^2"""
    # raiz = sqrt(2*(3 - a**2 * x**4))
    raiz = -6*xsqrta(x, a)*(2+a*x**2)+14*a*x**2 + 6
    return (2*x/3) + raiz/(6*a*x)

functions = {}
functions["bab"] = bab
functions["bab2"] = bab2
functions["stac"] = stac
functions["g11"] = g11
functions["g120"] = g120
functions["g121"] = g121
functions["g122"] = g122
functions["g123"] = g123
functions["g21"] = g21
functions["g220"] = g220
functions["g221"] = g221
functions["g222"] = g222
functions["g223"] = g223

npts = 12
a = mp.mpf(2)
x0 = 1/a
root = 1/sqrt(a)

vectors = {}
for i, key in enumerate(functions.keys()):
    vectors[key] = []
    vectors[key].append(mp.mpf(x0))


for j, (key, func) in enumerate(functions.items()):
    vector = vectors[key]
    for i in range(npts-1):
        try:
            fval = mp.mpf(func(vector[i], a))
            if mp.fabs(fval - root) > mp.fabs(vector[-1] - root):
                break
            vector.append(fval)
        except Exception as e:
            break
    for k in range(i):
        val = mp.mpf(vector[k]) - root
        val = mp.log10(mp.fabs(val))
        # val = mp.log10(mp.fabs(val))
        vector[k] = val
    plt.plot(vectors[key], label=key, marker=".")
    
plt.legend()
plt.show()