Random complex square-root function system

The random function system z -> ±√(z - omega), whose points settle on the Julia set of z² + omega.

• Name: inverse_julia
• Field: mathematics/fractals
• Formalism: discrete map
• References: Barnsley 1988; Dilão 2023

Select this system

from skcomplex.systems import create_simulation, create_system, describe

system = create_system('inverse_julia')
entry = describe('inverse_julia')
trajectory = create_simulation('inverse_julia').prepare().run().trajectory_

Factory

@register_system(
    field='mathematics/fractals',
    tags=('fractal', 'complex', 'iterated function system', 'random map', 'julia set'),
    references=('Barnsley 1988', 'Dilão 2023'),
    defaults={'until': 30, 'dt': 1},
)
def inverse_julia(
    omega: complex = -0.123 + 0.745j,
    *,
    p: float = 0.5,
    n_points: int = 5000,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the random function system `z -> ±√(z - omega)`, whose points settle on the Julia set of `z² + omega`.

    At every step each point moves to the positive square root of `z - omega` with probability `p` and to the
    negative one otherwise (Dilão 2023, Exercise 5.8). Both branches invert the quadratic map `z -> z² + omega`, which
    repels its Julia set, so the inverse branches attract the points onto the Julia set, drawn by the chaos game with
    weights `p` and `1 - p`. For `omega = 0` the Julia set is the unit circle, and the default is Douady's rabbit.

    Args:
        omega:
            The constant of the quadratic map.

        p:
            The probability of the positive square root.

        n_points:
            The number of points, started uniformly in the square of side 2 about the origin.

        random_state:
            The seed or generator of the initial points.

    Returns:
        system:
            The system, with the points `points.z`.

    Examples:
        >>> from skcomplex.systems.mathematics.fractals import inverse_julia
        >>> inverse_julia(0.0, n_points=10).elements[0].prepare().size_
        10
    """
    rng = np.random.default_rng(random_state)
    start = rng.uniform(-1.0, 1.0, n_points) + 1j * rng.uniform(-1.0, 1.0, n_points)
    values = {'z': start, 'omega': complex(omega), 'p': p}
    return System(elements=[GroupElement('points', element=InverseJuliaPointElement, size=n_points, values=values)])

Components

Elements:

• InverseJuliaPointElement, from Element

class InverseJuliaPointElement(Element):
    """Point of the inverse iteration of `z -> z² + omega`, which moves to a square root of `z - omega` drawn at random.

    The point moves to `+√(z - omega)` with probability `p` and to `-√(z - omega)` otherwise, with the principal
    square root.

    Args:
        label:
            The label of the element or the group.

        visual:
            How the element is drawn, or `None` for `default_visual_`.

        z:
            The initial point.

        omega:
            The constant of the quadratic map.

        p:
            The probability of the positive square root.

    Examples:
        >>> from skcomplex.systems.mathematics.fractals import InverseJuliaPointElement
        >>> InverseJuliaPointElement('point', p=0.3).p
        0.3
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        z: Annotated[object, StateSpec(domain=CC)] = 0.5,
        omega: Annotated[object, ParameterSpec(domain=CC)] = 0.0,
        p: Annotated[object, ParameterSpec(domain=RealSet.closed(0, 1))] = 0.5,
    ) -> None:
        super().__init__(label, visual=visual)
        self.z = z
        self.omega = omega
        self.p = p

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the move to the square root drawn in the step.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the element.
        """
        return [UpdateLaw(s.z, (2 * bernoulli('branch', s.p) - 1) * sqrt(s.z - s.omega))]

Equations

points.z[i].real: z_points_re(t + 1) == ((omega_points_im - z_points_im(t))^2 + (omega_points_re - z_points_re(t))^2)^(1/4)*(2*where(bernoulli_branch_points < p_points, 1, 0) - 1)*cos(1/2*arctan2(-omega_points_im + z_points_im(t), -omega_points_re + z_points_re(t)))
points.z[i].imag: z_points_im(t + 1) == ((omega_points_im - z_points_im(t))^2 + (omega_points_re - z_points_re(t))^2)^(1/4)*(2*where(bernoulli_branch_points < p_points, 1, 0) - 1)*sin(1/2*arctan2(-omega_points_im + z_points_im(t), -omega_points_re + z_points_re(t)))