Skip to content

polykin.math¤

roots_xtanx cached ¤

roots_xtanx(
    a: float, N: int, xtol: float = 1e-06
) -> FloatVector

Determine the first N roots of the transcendental equation x_n*tan(x_n) = a.

PARAMETER DESCRIPTION
a

Parameter.

TYPE: float

N

Number of roots.

TYPE: int

xtol

Tolerance.

TYPE: float DEFAULT: 1e-06

RETURNS DESCRIPTION
FloatVector

Roots of the equation.

Examples:

Determine the first 4 roots for a=1.

>>> from polykin.math import roots_xtanx
>>> roots_xtanx(1.0, 4)
array([0.86033359, 3.42561846, 6.43729818, 9.52933441])
Source code in src/polykin/math/special.py 
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
@functools.cache
def roots_xtanx(a: float, N: int, xtol: float = 1e-6) -> FloatVector:
    r"""Determine the first `N` roots of the transcendental equation 
    `x_n*tan(x_n) = a`.

    Parameters
    ----------
    a : float
        Parameter.
    N : int
        Number of roots.
    xtol : float
        Tolerance.

    Returns
    -------
    FloatVector
        Roots of the equation.

    Examples
    --------
    Determine the first 4 roots for a=1.
    >>> from polykin.math import roots_xtanx
    >>> roots_xtanx(1.0, 4)
    array([0.86033359, 3.42561846, 6.43729818, 9.52933441])
    """
    result = np.zeros(N)

    if a == 0:
        for i in range(N):
            result[i] = i*pi
    elif a == inf:
        for i in range(N):
            result[i] = (i + 1/2)*pi
    elif a < 1e2:
        for i in range(N):
            sol = fzero_newton(f=lambda x: x*tan(x) - a,
                               x0=pi*(i + 0.5*a/(a+1)),
                               xtol=xtol,
                               ftol=1e-10)
            result[i] = sol.x
    else:
        for i in range(N):
            x0 = i*pi
            e = 1.0
            for _ in range(0, 10):
                e_new = arctan(a/(x0 + e))
                if abs(e_new - e) < xtol:
                    break
                e = e_new
            result[i] = x0 + e

    return result