Interpolation and Fixed Times

Mix.install([
  {:integrator, github: "woodward/integrator"},
  {:kino_vega_lite, "~> 0.1"}
])

Interpolation

By default, refine: 4 for integrator DormandPrince45. This means that four points are interpolated for every solution of the ODE.

alias Integrator.SampleEqns
alias Integrator.DataCollector

t_initial = Nx.f64(0.0)
t_final = Nx.f64(20.0)
x_initial = Nx.f64([2.0, 0.0])
{:ok, pid} = DataCollector.start_link()
output_fn = &DataCollector.add_data(pid, &1)

opts = [type: :f64, output_fn: output_fn]

_solution = Integrator.integrate(&SampleEqns.van_der_pol_fn/2, t_initial, t_final, x_initial, opts)
points = DataCollector.get_data(pid)

Visualizing refine: 4:

alias VegaLite, as: VLq
alias Integrator.Point

defmodule VanDerPol do
  def plot(points) do
    data =
      points
      |> Enum.map(&Point.to_number(&1))
      |> Enum.map(fn point ->
        %{t: t, x: x} = point

        [
          %{t: t, x: List.first(x), x_value: "x[0]"},
          %{t: t, x: List.last(x), x_value: "x[1]"}
        ]
      end)
      |> List.flatten()

    VL.new(
      width: 600,
      height: 400,
      title: "Solution of van der Pol Equation (μ = 1) with Dormand-Prince45"
    )
    |> VL.mark(:line, point: true, tooltip: true)
    |> VL.encode_field(:x, "t", type: :quantitative)
    |> VL.encode_field(:y, "x", type: :quantitative)
    |> VL.encode_field(:color, "x_value", type: :nominal)
    |> VL.data_from_values(data)
    |> Kino.VegaLite.new()
    |> Kino.render()
  end
end

VanDerPol.plot(points)

You can turn off interpolation by setting refine: 1:

t_initial = Nx.f64(0.0)
t_final = Nx.f64(20.0)
x_initial = Nx.f64([2.0, 0.0])

{:ok, pid} = DataCollector.start_link()
output_fn = &DataCollector.add_data(pid, &1)

opts = [type: :f64, output_fn: output_fn, refine: 1]


_solution =
  Integrator.integrate(&SampleEqns.van_der_pol_fn/2, t_initial, t_final, x_initial, opts)

points = DataCollector.get_data(pid)
VanDerPol.plot(points)

Note how much “chunkier” the plot is without the interpolated points. These points are solely those from the Runge-Kutta simulation. Note that these values can also be accessed by solution.ode_t and solution.ode_x.

Finally, you can output data at fixed times. For example, let’s print out data at 0.1 second intervals:

t_initial = Nx.f64(0.0)
t_final = Nx.f64(20.0)
x_initial = Nx.f64([2.0, 0.0])

{:ok, pid} = DataCollector.start_link()
output_fn = &DataCollector.add_data(pid, &1)

opts = [type: :f64, refine: 1, fixed_output_times?: Nx.u8(1), fixed_output_step: Nx.f64(0.1), output_fn: output_fn]

_solution = Integrator.integrate(&SampleEqns.van_der_pol_fn/2, t_initial, t_final, x_initial, opts)
points = DataCollector.get_data(pid)
VanDerPol.plot(points)