specular.ode.euler_scheme_1

Copy

euler_scheme_1(F, t_0, T, u_0, u_1, *, n_steps)

Apply the explicit two-step specular Euler method of Type 1 to a scalar initial-value problem.

Parameters

Name	Type	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.
u_1	real scalar	Finite, externally computed starting value at the first time node t[1].
n_steps	positive integer, keyword-only	Number of uniform time intervals.

Returns

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

Notes

With \(h=(T-t_0)/N\), \(N=\texttt{n_steps}\), and \(\mathcal C=\mathcal C_1\), the recurrence is

\[ u_{n+1}=u_n+h\mathcal C\!\left( F(t_n,u_n),F(t_{n-1},u_{n-1})\right),\qquad n=1,\ldots,N-1. \]

The solver keeps the supplied u_1 unchanged. For n_steps=1, it returns u_0 and u_1 without evaluating F. The all-ones sigma array describes the unscaled-mean convention, including the externally supplied first interval; it does not describe how the caller produced u_1.

This method is generically first-order consistent. A more accurate starter does not give an unconditional higher-order theorem. See Type 1. Angular means use the selected backend; the outer iteration remains scalar and the stored arrays use float64.

Examples

Copy

import specular.ode

def F(t, u):
    return -u

t_0, T, u_0 = 0.0, 1.0, 1.0
n_steps = 100
h = (T - t_0) / n_steps
u_1 = u_0 + h * F(t_0, u_0)
result = specular.ode.euler_scheme_1(
    F, t_0, T, u_0, u_1, n_steps=n_steps,
)
print(result.u[-1])

Source code in specular/ode/_euler.py

Copy

def euler_scheme_1(
    F: ScalarField,
    t_0: RealScalar,
    T: RealScalar,
    u_0: RealScalar,
    u_1: RealScalar,
    *,
    n_steps: int,
) -> ODEResult:
    r"""Apply the explicit two-step specular Euler method of Type 1.

    ``u_1`` is an externally supplied value at the first represented node.
    """

    t_values, step_sizes, initial_value, first_value = _two_step_inputs(
        F,
        t_0,
        T,
        u_0,
        u_1,
        n_steps,
    )
    counted_field = _FieldEvaluationCounter(F)
    step_count = len(step_sizes)
    u_values = np.empty(step_count + 1, dtype=np.float64)
    u_values[0] = initial_value
    u_values[1] = first_value

    if step_count > 1:
        F_previous = _field_value(
            counted_field,
            float(t_values[0]),
            initial_value,
            step=0,
        )
        for step in range(1, step_count):
            u_n = float(u_values[step])
            F_current = _field_value(
                counted_field,
                float(t_values[step]),
                u_n,
                step=step,
            )
            slope = float(scaled_mean(F_current, F_previous))
            u_values[step + 1] = _advance(
                u_n,
                float(step_sizes[step]),
                slope,
                step=step,
            )
            F_previous = F_current

    return ODEResult(
        t=t_values,
        u=u_values,
        sigma=np.ones(step_count, dtype=np.float64),
        number_of_field_evaluations=counted_field.number_of_field_evaluations,
    )