defmodule Fluxion.ODE.Euler do
  @moduledoc """
  The Euler method for solving the initial value problem for ordinary
  differential equations of the form @y'(t) = f(t, y(t))@.

  Euler's method defines a series of @y@ values according to the following:

  @y\\_{k+1} = y\\_k + hf(t\\_k, y\\_k)@ with @k = 0, 1, ..., n@

  In this module, the value for @y\\_0@ will be referred to as "initial" and
  the value for @h@ will be referred to as the "stride".

  The stride is a crucial parameter for Euler's method. Euler's method is
  convergent because as @h \\to 0@ the series of values @y\\_k@ approaches
  @y(t\\_k)@.

  Convergence rate is @O(h)@ and numeric error accumulates as @O(1/h)@,
  so properly selected stride is important to balance convergence and
  numeric error.
  """

  @doc """
  Step Euler's solution @n@ times.

  ## Example

  For the purpose of example, start from a known function for @y@ like
  @y(t) = sin(t)@.
  This gives @y'(t) = cos(t)@.
  Given this, it's possible to approximate the real solution with
  Euler's method.

      iex> defmodule TSq do
      ...>   @behaviour Fluxion.ODE
      ...>   @impl true
      ...>   def f(x, _y_k), do: :math.cos(x)
      ...> end
      iex> steps = 10
      iex> stride = 0.01
      iex> # Compute the real value of y(t) after 10 steps
      iex> real_y_1 = :math.sin(stride*steps)
      iex> # Approximate the value of y(t) after 10 steps
      iex> approx_y_1 = Fluxion.ODE.Euler.step(TSq, 0 , stride, steps)
      iex> # Compare the results
      iex> Fluxion.Math.is_equal(real_y_1, approx_y_1, 1.0e-3)
      :true
  """
  @spec step(
          ode :: Fluxion.ODE.t(),
          initial :: number,
          stride :: number,
          n :: number
        ) :: number
  def step(ode, initial, stride, n \\ 1) do
    iterate(ode, initial, stride, 0, n)
  end

  # Handle the last iteration
  defp iterate(ode, y_k, stride, t, 1) do
    y_k + stride * Fluxion.ODE.f!(ode, t, y_k)
  end

  # Handle all other iteration steps
  defp iterate(ode, y_k, stride, t, n) do
    dydt = Fluxion.ODE.f!(ode, t, y_k)
    y_next = y_k + stride * dydt
    iterate(ode, y_next, stride, t + stride, n - 1)
  end
end