specular.ode.ellipse_scheme

Copy

ellipse_scheme(F, t_0, T, u_0, *, n_steps, sigma_n=None, third_order=False, fourth_order=False, minimize_defect=False, derivatives_of_F=None, derivative_step=None, atol=1e-12, rtol=1e-10, max_iter=100)

Solve a scalar initial-value problem on [t_0, T] with the implicit specular ellipse method.

Parameters

Name	Type and default	Description
F	callable	Field F(t, u) returning a finite real scalar.
t_0	real scalar	Finite initial time.
T	real scalar	Finite final time, strictly greater than t_0.
u_0	real scalar	Finite initial value at t_0.
n_steps	positive integer, keyword-only	Number of time intervals; the result has n_steps + 1 nodes.
sigma_n	positive real scalar or callable, default None	Required in the base mode. A callable sigma_n(n, t_n, u_n, h) returns a finite positive scale and is evaluated once at the accepted left endpoint. Here h is the represented interval t[n + 1] - t[n]. Must be None in an automatic mode.
third_order	bool, default False	Select a positive scale using numerical left-endpoint defect cancellation.
fourth_order	bool, default False	Use the coupled two-endpoint scale in cases E5a/E5b and scale 1.0 in the other cases.
minimize_defect	bool, default False	Apply the full E1--E6 two-endpoint defect-minimization rule, including zero- and infinite-scale limits.
derivatives_of_F	callable or None, default None	In an automatic mode, optionally map the NumPy array [t, u] to [L_F F, L_F^2 F], also of shape (2,). Omit to estimate these derivatives numerically.
derivative_step	positive real scalar or None, default None	Finite-difference step for automatic modes. Cannot be combined with derivatives_of_F.
atol	real scalar, default 1e-12	Finite nonnegative absolute tolerance for the nonlinear solve.
rtol	real scalar, default 1e-10	Finite nonnegative relative tolerance for the nonlinear solve. At least one of atol and rtol must be positive.
max_iter	positive integer, default 100	Limit for each local nonlinear iteration or root-refinement phase. Bracket searches may make additional field evaluations.

All parameters after u_0 are keyword-only. The three automatic-mode flags are mutually exclusive. Derivative controls are valid only in an automatic mode. Here \(L_F=\partial_t+F\partial_u\) denotes differentiation along solution curves.

Returns

ODEResult: NumPy float64 arrays t and u of shape (n_steps + 1,), a scale array sigma of shape (n_steps,), and the number of solver calls to F.

Notes

For a prescribed scale, the update is

\[ u_{n+1}=u_n+h\mathcal C_{\sigma_n}\!\left( F(t_{n+1},u_{n+1}),F(t_n,u_n)\right). \]

The scale is frozen during each implicit solve. Angular means use the selected calculation backend; the outer iteration uses Python floats and NumPy float64 arrays. See Prescribed scale.

Local nonlinear convergence is conditional. A failed solve raises RuntimeError; a smaller time step may help. Exceptions from F are preserved. See Nonlinear solve controls.

The order flags do not guarantee their named order for arbitrary fields. Third-order convergence requires a sufficiently smooth, bounded selected branch and accurate derivatives. The fourth-order theorem additionally requires the uniform E5a-or-E5b condition near the solution and the stated nondegeneracy assumptions; steps using the scale-1 fallback do not inherit that guarantee. minimize_defect is a defect-minimization mode with no fourth-order or maximal-order guarantee. See Automatic scale-selection modes.

Numerical derivative estimation may sample F outside [t_0, T]; the field must be defined on a neighborhood of the interval, or exact derivatives must be supplied. Only minimize_defect records the result sentinels 0.0 and inf for limiting methods. They are not valid prescribed scales or public scaled_mean scale arguments.

Examples

Copy

import specular.ode

result = specular.ode.ellipse_scheme(
    lambda t, u: -u,
    0.0,
    1.0,
    1.0,
    n_steps=100,
    sigma_n=1.0,
)
print(result.t[-1], result.u[-1])

Source code in specular/ode/_ellipse.py

Copy

def ellipse_scheme(
    F: ScalarField,
    t_0: RealScalar,
    T: RealScalar,
    u_0: RealScalar,
    *,
    n_steps: int,
    sigma_n: RealScalar | StepScale | None = None,
    third_order: bool = False,
    fourth_order: bool = False,
    minimize_defect: bool = False,
    derivatives_of_F: VectorToVectorFunc | None = None,
    derivative_step: RealScalar | None = None,
    atol: RealScalar = 1e-12,
    rtol: RealScalar = 1e-10,
    max_iter: int = 100,
) -> ODEResult:
    r"""Apply the scalar specular ellipse method on ``[t_0, T]``.

    In the base mode, ``sigma_n`` is a positive scalar or a callable
    ``sigma_n(n, t_n, u_n, h)`` and is frozen during its implicit step. The
    fourth argument is the represented step size ``t[n + 1] - t[n]``.
    ``third_order=True`` numerically enforces left-endpoint defect
    cancellation. ``fourth_order=True`` couples the endpoint update to the
    two-endpoint scale in cases E5a and E5b and otherwise uses scale ``1.0``.
    ``minimize_defect=True`` instead applies the full E1--E6 two-endpoint
    defect-minimization rule, including its zero- and infinite-scale limits.

    When ``derivatives_of_F`` is omitted, ``L_F F`` and ``L_F^2 F`` are
    estimated from ``F`` along a local RK4 flow. Otherwise the callback must
    map ``[t, u]`` to ``[L_F F, L_F^2 F]``.
    """

    (
        initial_time,
        final_time,
        initial_value,
        step_count,
        selected_sigma,
        use_third,
        use_fourth,
        use_minimize,
        selected_derivatives,
        selected_derivative_step,
        absolute_tolerance,
        relative_tolerance,
        iteration_limit,
    ) = _validated_inputs(
        F,
        t_0,
        T,
        u_0,
        n_steps=n_steps,
        sigma_n=sigma_n,
        third_order=third_order,
        fourth_order=fourth_order,
        minimize_defect=minimize_defect,
        derivatives_of_F=derivatives_of_F,
        derivative_step=derivative_step,
        atol=atol,
        rtol=rtol,
        max_iter=max_iter,
    )
    counted_field = _FieldEvaluationCounter(F)

    t_values, step_sizes = _time_grid(
        initial_time,
        final_time,
        step_count,
    )
    if (
        selected_derivative_step is not None
        and selected_derivatives is None
    ):
        derivative_times = (
            t_values
            if use_fourth or use_minimize
            else t_values[:-1]
        )
        for sample_time_value in derivative_times:
            sample_time = float(sample_time_value)
            for radius in (
                0.5 * selected_derivative_step,
                selected_derivative_step,
            ):
                forward_time = sample_time + radius
                backward_time = sample_time - radius
                if (
                    not math.isfinite(forward_time)
                    or not math.isfinite(backward_time)
                    or forward_time == sample_time
                    or backward_time == sample_time
                ):
                    raise ValueError(
                        "derivative_step does not produce distinct float64 "
                        "samples on [t_0, T]"
                    )

    u_values = np.empty(step_count + 1, dtype=np.float64)
    sigma_values = np.empty(step_count, dtype=np.float64)
    u_values[0] = initial_value
    previous_sigma: float | None = None

    for step in range(step_count):
        t_n = float(t_values[step])
        t_next = float(t_values[step + 1])
        u_n = float(u_values[step])
        h_n = float(step_sizes[step])

        if not use_third and not use_fourth and not use_minimize:
            if callable(selected_sigma):
                sigma = _positive_scale(
                    selected_sigma(step, t_n, u_n, h_n), step=step
                )
            else:
                assert selected_sigma is not None
                sigma = selected_sigma
            F_left = _field_value(counted_field, t_n, u_n, step=step)
            u_next = _fixed_scale_step(
                counted_field,
                step=step,
                t_next=t_next,
                u_n=u_n,
                h_n=h_n,
                sigma=sigma,
                F_left=F_left,
                atol=absolute_tolerance,
                rtol=relative_tolerance,
                max_iter=iteration_limit,
            )
        else:
            left = _field_quantities(
                counted_field,
                t_n,
                u_n,
                selected_derivatives,
                selected_derivative_step,
                step=step,
            )
            if use_third:
                sigma = _third_order_scale(
                    left, previous_sigma, step=step
                )
                u_next = _fixed_scale_step(
                    counted_field,
                    step=step,
                    t_next=t_next,
                    u_n=u_n,
                    h_n=h_n,
                    sigma=sigma,
                    F_left=left[0],
                    atol=absolute_tolerance,
                    rtol=relative_tolerance,
                    max_iter=iteration_limit,
                )
            else:
                u_next, sigma = _coupled_step(
                    counted_field,
                    step=step,
                    t_next=t_next,
                    u_n=u_n,
                    h_n=h_n,
                    left=left,
                    derivatives_of_F=selected_derivatives,
                    derivative_step=selected_derivative_step,
                    minimize_defect=use_minimize,
                    atol=absolute_tolerance,
                    rtol=relative_tolerance,
                    max_iter=iteration_limit,
                )

        u_values[step + 1] = u_next
        sigma_values[step] = sigma
        if use_third:
            previous_sigma = sigma

    return ODEResult(
        t=t_values,
        u=u_values,
        sigma=sigma_values,
        number_of_field_evaluations=counted_field.number_of_field_evaluations,
    )