ex7: problem 8

exercise7/main.typ

@@ -41,3 +41,7 @@ this document was created using
 = problem 4
 
 #include "problem4.typ"
+
+= problem 8
+
+#include "problem8.typ"

exercise7/problem8.typ

@@ -0,0 +1,164 @@
+#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.