Harmonic Oscillator

Let's demonstrate basic capabilities of SciLean on a simple program simulating harmonic oscillator. The energy/Hamiltonian of harmonic oscillator is

def H (m k : Float) (x p : Float) := (1/(2*m)) * p^2 + k/2 * x^2

Once we have an energy of a physical system we can derive its dynamics using Hamilton's equations

$$\begin{align} \dot x &= \frac{\partial H}{\partial p} \\ \dot p &= -\frac{\partial H}{\partial x} \end{align} $$

In SciLean, we define function solver that approximatelly solves these equations with the following command

approx solver (m k : Float)
  := odeSolve (fun (t : Float) (x,p) => ( ∇ (p':=p), H m k x  p',
                                         -∇ (x':=x), H m k x' p))
by
  ...

The expression after := is just a specification of what the function solver is supposed to approximate and the actual definition is given where the dots ... are.

The general idea behind SciLean is that we first specify what we want to compute and only then we describe how we do it. In this case, the function SciLean.odeSolve is non-computable function that can't be executed but has well defined meaning of returning the solution of given differential equation. In the place of the dots ... we will manipulate this function symbolically and transform the whole specification into an expression that can be evaluated approximatelly.

Let's have a look how we can turn the specification into implementation

open SciLean
approx solver (m k : Float)
  := odeSolve (fun (t : Float) (x,p) => ( ∇ (p':=p), H m k x  p',
                                         -∇ (x':=x), H m k x' p))
by
  unfold Hm : Float
k : Float
⊢ Approx (Filter.atTop.prod (?m.11758 m k))
  (odeSolve fun t x =>
    match x with
    | (x, p) => (∇ (p':=p), (1 / (2 * m) * p' ^ 2 + k / 2 * x ^ 2), -∇ (x':=x), (1 / (2 * m) * p ^ 2 + k / 2 * x' ^ 2)))  -- unfold definition of H
  autodiffm : Float
k : Float
⊢ let ydf := (2 * m)⁻¹;
let ydf_1 := k / 2;
Approx (Filter.atTop.prod (?m.11758 m k))
  (odeSolve fun t x =>
    let ydf_2 := x.2;
    let ydf_3 := x.1;
    (ydf * (2 * ydf_2), -(ydf_1 * (2 * ydf_3))))  -- run automatic differentiation

  -- apply RK4 method
  simp_rw (config:={zeta:=false}) [odeSolve_fixed_dt rungeKutta4 sorry_proofm : Float
k : Float
⊢ let ydf := (2 * m)⁻¹;
let ydf_1 := k / 2;
Approx (?m.3965 m k) fun t₀ t x₀ =>
  limit n → ∞,
    odeSolveFixedStep
      (rungeKutta4 fun t x =>
        let ydf_2 := x.2;
        let ydf_3 := x.1;
        (ydf * (2 * ydf_2), -(ydf_1 * (2 * ydf_3))))
      n t₀ t x₀]

  -- approximate the limit with a fixed `n`
  approx_limit n sorry_proofm : Float
k : Float
n : ℕ
⊢ Approx (?m.11758 m k) fun t₀ t x₀ =>
  odeSolveFixedStep
    (rungeKutta4 fun t x =>
      let ydf := x.2;
      let ydf_1 := x.1;
      ((2 * m)⁻¹ * (2 * ydf), -(k / 2 * (2 * ydf_1))))
    n t₀ t x₀
⊢ Float → Float → Type
⊢ (m k : Float) → Filter (?m.11756 m k)

There are two main steps in turning the specification into implementation. We need to compute the derivatives ∇ (p':=p), H m k x p' and ∇ (x':=x), H m k x' p and we need to approximate the function SciLean.odeSolve using some time stepping scheme.

To compute the derivatives, we unfold the definition of the Hamiltonian H and call autodiff tactic.

To approximate SciLean.odeSolve we rewrite the specification using the theorem SciLean.odeSolve_fixed_dt saying that SciLean.odeSolve can be approximated by a time stepping scheme, we choose Runge-Kutta 4 SciLean.rungeKutta4. j The theorem SciLean.odeSolve_fixed_dt produces a specification of the form limit n → ∞, .... We can decide to approximate this limit by taking a concrete value n and hope it is sufficiently large such that the approximation is good enough. To do this we call the tactic approx_limit and we end up with a goal that fully computable expression.

Now we have a function solver that takes mass m, stiffness k, number of substeps substep, initial time t₀, final time t, initial position x and momentum p and produces new position and momentum at time t.

fun m k substep t₀ t x =>
  match x with
  | (x, p) =>
    ApproxSolution.val (solver m k) substep t₀ t
      (x, p) : Float → Float → ℕ × Unit → Float → Float → Float × Float → Float × Float#check fun m k substep t₀ t (x,p) => solver m k substep t₀ t (x,p)

If we check solver directly

solver (m k : Float) :
  Approx (Filter.atTop.prod ⊤)
    (odeSolve fun t x =>
      match x with
      | (x, p) => (∇ (p':=p), H m k x p', -∇ (x':=x), H m k x' p))#check solver

we see that its type contains @Approx ?l ?spec where ?l specifies what are the parameters controling the approximation and ?spec is the specification of what is solver approximating.

Here is an example how to run solver

def sim : IO Unit := do

  let substeps := 1
  let m := 1.0
  let k := 10.0

  let Δt := 0.1
  let x₀ := 1.0
  let p₀ := 0.0
  let mut t := 0
  let mut (x,p) := (x₀, p₀)

  for _ in [0:50] do

    (x, p) := solver m k substeps t (t+Δt) (x, p)
    t += Δt

    -- poor man's visualization
    for (j : Nat) in [0:20] do
      if j.toFloat < 10*(x+1.0) then
        IO.print "o"
    IO.println ""

When we execute it we obtain expected sinusoidal wave

oooooooooooooooooooo
ooooooooooooooooooo
oooooooooooooooo
oooooooooooooo
oooooooooo
ooooooo
ooooo
oo
o
o
o
ooo
ooooo
oooooooo
ooooooooooo
oooooooooooooo
ooooooooooooooooo
ooooooooooooooooooo
oooooooooooooooooooo
oooooooooooooooooooo
oooooooooooooooooooo
oooooooooooooooooo
oooooooooooooooo
ooooooooooooo
oooooooooo
ooooooo
oooo
oo
o
o
o
ooo
ooooo
oooooooo
ooooooooooo
oooooooooooooo
ooooooooooooooooo
ooooooooooooooooooo
oooooooooooooooooooo
oooooooooooooooooooo
oooooooooooooooooooo
oooooooooooooooooo
oooooooooooooooo
ooooooooooooo
oooooooooo
ooooooo
oooo
oo
o
o
#eval sim