165 lines
3.5 KiB
Typst
165 lines
3.5 KiB
Typst
#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.
|