#import "@preview/physica:0.9.6": *
#import "@preview/unify:0.7.1": qty
#import "@preview/cetz:0.3.2"
#import "@preview/cetz-plot:0.1.1"

== a)

we are given a system of ODEs describing the movement of a pendulum
$
  cases(
    theta''(t) = - g/L sin(theta(t)),
    theta(0) = pi/4,
    theta'(0) = 0
  )
$
where $g = qty("9.81", "m/s^2")$ and $L = qty("1", "m")$.

let $theta_k = theta^((k))$ in
$
  cases(
    theta_0(0) = pi/4,
    theta_1(0) = 0,
    theta_2(t) = - g/L sin(theta_0(t)),
    theta'_0 = theta_1,
    theta'_1 = theta_2
  )
$


== b)

recall euler's method: $y_(n+1) = y_n + h y'_n$

#let code = ```typc
  let h = 0.01; let g = 9.81; let L = 1;

  let t = 0;
  let theta_0 = calc.pi / 4;
  let theta_1 = 0;
  let data = ();
  while t < 5 {
    data.push((t, theta_0, theta_1));
    theta_0 += h * theta_1;
    theta_1 += h * (-g/L * calc.sin(theta_0));
    t += h;
  };
  data
```

#let euler-data = eval(code.text.split("\n").join())

$theta$ is at $t = 1$ approximately #calc.round(euler-data.at(100).at(1), digits: 4) (with step size $h = 0.01$ so i can reuse the code later).

#code


== c)

recall heun's method:
$
  y_(n+1) = y_n + h/2 [f(t_n, y_n) + f(t_(n+1), y_n + h f(t_n, y_n))]
$

#let code = ```typc
  let h = 0.01; let g = 9.81; let L = 1;

  let t = 0;
  let theta_0 = calc.pi / 4;
  let theta_1 = 0;
  let data = ();
  while t < 5 {
    data.push((t, theta_0, theta_1));
    let euler_theta_1 = h * (-g/L
            * calc.sin(theta_0 + h * theta_1));
    theta_0 += h/2 * (theta_1 + euler_theta_1);
    theta_1 += h * (-g/L * calc.sin(theta_0));
    t += h;
  };
  data
```

#let heun-data = eval(
  code
    .text
    .split("\n")
    .filter(l => {
      if l.len() > 2 {
        l.rev().trim().rev().slice(0, 2) != "//"
      } else { true }
    })
    .join(),
)

$theta$ is approxmately equal to #calc.round(heun-data.at(100).at(1), digits: 4) at $t = 1$ (with step size $h = 0.01$, so i can reuse the code later).

#code


== d)


#cetz.canvas({
  import cetz.draw: *
  import cetz-plot: *
  plot.plot(
    axis-style: "school-book",
    x-label: $t$,
    y-label: $theta$,
    size: (10, 8),
    x-tick-step: 0.5,
    y-tick-step: 0.5,
    {
      plot.add(euler-data, label: "euler")
      plot.add(heun-data, label: "heun")
    },
  )
})

I'm thought heun's method was wrongly implemented, because I didn't bother
using explicit function definitions, which complicates things, but it seems to
be more correct than euler's method.


== e)

using the function $E(theta_0, theta_1) = 1/2 L^2 (theta_1)^2 + g L (1
  - cos(theta_0))$ I can calculate the energy levels for both graphs at each
point in time.

#{
  let g = 9.81
  let L = 1
  let E(theta_0, theta_1) = (
    1 / 2 * L * L * theta_1 * theta_1 + g * L * (1 - calc.cos(theta_0))
  )
  let heun-energy = heun-data.map(data => (
    data.at(0),
    E(data.at(1), data.at(2)),
  ))
  let euler-energy = euler-data.map(data => (
    data.at(0),
    E(data.at(1), data.at(2)),
  ))
  cetz.canvas({
    import cetz.draw: *
    import cetz-plot: *
    plot.plot(
      axis-style: "school-book",
      x-label: $t$,
      y-label: $theta$,
      size: (10, 8),
      x-tick-step: 0.5,
      y-tick-step: 0.5,
      {
        plot.add(euler-energy, label: "euler")
        plot.add(heun-energy, label: "heun")
      },
    )
  })
}

we can see that neither are constant. this is bad (probably). however, euler's
method seems to be most correct in this sense too, if we are modelling a system
devoid of friction, the energy should remain the same.