Damped oscillations

This example is borrowed from owl-ode.

open Owl_ode
open Owl_ode.Types
open Owl_plplot

let damped_noforcing a (xs, ps) _ : Owl.Mat.mat =
  Owl.Mat.((xs *$ -1.0) + (ps *$ (-1.0 *. a)))

let a = 1.0
let dt = 0.1

let plot_sol fname t sol1 sol2 sol3 =
  let open Owl in
  let h = Plot.create fname in
  let open Plot in
  set_foreground_color h 0 0 0;
  set_background_color h 255 255 255;
  set_title h "Damped oscillations" ;
  plot ~h ~spec:[ RGB (0, 0, 255); LineStyle 1 ] t (Mat.col sol1 0);
  plot ~h ~spec:[ RGB (0, 255, 0); LineStyle 1 ] t (Mat.col sol2 0);
  plot ~h ~spec:[ RGB (255, 0, 0); LineStyle 1 ] t (Mat.col sol3 0);
  (* XXX: I could not figure out how to make the legend black instead of red *)
  legend_on h ~position:NorthEast [| "Leapfrog"; "Ruth3"; "Symplectic Euler" |];
  output h
;;
val damped_noforcing :
  float ->
  (float, Bigarray.float64_elt) Owl_dense_matrix_generic.t *
  (float, Bigarray.float64_elt) Owl_dense_matrix_generic.t ->
  'a -> Owl.Mat.mat = <fun>
val a : float = 1.
val dt : float = 0.1
val plot_sol :
  string ->
  Owl_plplot.Plot.dsmat -> Owl.Mat.mat -> Owl.Mat.mat -> Owl.Mat.mat -> unit =
  <fun>

let _ =
  Scirep.Show.png (fun fn ->
    let x0 = Owl.Mat.of_array [| -0.25 |] 1 1 in
    let p0 = Owl.Mat.of_array [| 0.75 |] 1 1 in
    let t0, duration = 0.0, 15.0 in
    let f = damped_noforcing a in
    let tspec = T1 { t0; duration; dt } in
    let t, sol1, _ = 
      Ode.odeint (module Symplectic.D.Leapfrog) f (x0, p0) tspec () in
    let _, sol2, _ = 
      Ode.odeint Symplectic.D.ruth3 f (x0, p0) tspec () in
    let _, sol3, _ =
      Ode.odeint (module Symplectic.D.Symplectic_Euler) f (x0, p0) tspec ()
    in
    (* XXX: I'd prefer t to be already an Owl array as well *)
    plot_sol fn t sol1 sol2 sol3
)
;;

- : unit = ()