diff --git a/pendulum.md b/pendulum.md new file mode 100644 index 0000000..de205dc --- /dev/null +++ b/pendulum.md @@ -0,0 +1,225 @@ +--- +jupyter: + jupytext: + formats: ipynb,md + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.18.1 + kernelspec: + display_name: Julia 1.12 + language: julia + name: julia-1.12 +--- + +# Pendulum + +## Equations of motion +The dimensional equation for a pendulum with the constraints +- no friction +- rigid rod +- vertically movable anchor point with prescribed acceleration $a_0(t)$ with amplitude $c^{*}$ + +reads +$$ +\frac{\text{d}^2\theta}{\text{d} t^{*2}} = - \frac{a_0(t^{*})}{l^{*}} \sin(\theta) - \frac{g^{*}}{l^{*}} \sin(\theta) +$$ + +In the following, we will solely be concerned with a prescribed sinusoidal motion. Its acceleration is given by +$$ +a_0(t^{*}) = c^{*} \sin(k^{*} t^{*}) +$$ +with amplitude $c^{*}$ and wave number $k^{*}$. + +For non-dimensionalisation, we introduce a time scale $T$, which yields +$$ +\frac{1}{T^2}\frac{\text{d}^2\theta}{\text{d} t^2} = - \frac{c^{*}\sin(k^{*} T t) + g^{*}}{l^{*}} \sin(\theta). +$$ +Moreover, the following non-dimensional parameters are introduced. +$$ \tau := T^{-1} \sqrt{l^{*}/ g^{*}} \quad\text{(gravitational time scale)}$$ +$$ c := \frac{c^{*}}{g^{*}} \quad\text{(relative anchor point acceleration)}$$ + +$$ k := k^{*}T \quad\text{(anchor point oscillation time scale)}$$ + +The non-dimensional equation then reads +$$ +\frac{\text{d}^2\theta}{\text{d} t^{2}} = - \frac{1}{\tau^2} (c \sin(k t) + 1) \sin(\theta), +$$ +which is consistent with the Buckingham theorem given five dimensional variables +($t^*$,$g^*$,$c^*$,$k^*$,$l^*$) and two dimensions (time, length). + +## Numerical method +We define the angular velocity of the pendulum $\Omega$ to reformulate our problem into a set of two first-order ordinary +differential equations, +$$ +\frac{\text{d}\Omega}{\text{d}t} = - \frac{1}{\tau^2} (c \sin(k t) + 1) \sin(\theta), +$$ +$$ +\text{d}\theta / \text{d}t = \Omega, +$$ +which can be solved numerically. + +We employ the classical Runge-Kutta method (RK4) for $\text{d}\mathbf{x}/\text{d}t = f(\mathbf{x},t)$, i.e. +$$ \mathbf{x}_{n+1} = \mathbf{x}_n + \frac{\Delta t}{6} ( r_1 + 2 r_2 + 2 r_3 + r_4 ),$$ +$$ r_1 = f(t_n, \mathbf{x}_n),$$ +$$ r_2 = f(t_n + \Delta t /2, \mathbf{x}_n + r_1 \Delta t /2),$$ +$$ r_3 = f(t_n + \Delta t /2, \mathbf{x}_n + r_2 \Delta t /2),$$ +$$ r_4 = f(t_n + \Delta t, \mathbf{x}_n + r_3 \Delta t),$$ +where $\mathbf{x} := (\theta,\Omega,t)^T$ is the state vector of our problem. +This enables us to get a discrete-in-time map of the form +$$ \mathbf{x}_{n+1} = f(\mathbf{x}_n,\mathbf{p}) $$ +with $\mathbf{x}$ and $\mathbf{p}$ being the state and parameter vectors, respectively. Note that $\mathbf{p}$ contains +the parameters of the physical system, as well as the numerical parameter ($\Delta t$) since both affect the evolution of +the discrete numerical system. + +```julia +struct State + θ :: Float64 + Ω :: Float64 + t :: Float64 +end + +struct Parameter + τ :: Float64 + c :: Float64 + k :: Float64 + Δt :: Float64 +end + +function dsmap(x::State,p::Parameter) + f = (x::State,p::Parameter) -> -(1.0/p.τ^2)*(p.c*sin(p.k*x.t)+1)*sin(x.θ) + coeff = [p.Δt/6.0,p.Δt/3.0,p.Δt/3.0,p.Δt/6.0] + step = [0.0,p.Δt/2.0,p.Δt/2.0,p.Δt] + rΩ = f(x,p) + rθ = x.Ω + Ωnp1 = x.Ω + coeff[1]*rΩ + θnp1 = x.θ + coeff[1]*rθ + for ii = 2:length(coeff) + rΩ = f(State(x.θ+rθ*step[ii],NaN,x.t+step[ii]),p) + rθ = x.Ω + rΩ*step[ii] + θnp1 += coeff[ii]*rθ + Ωnp1 += coeff[ii]*rΩ + end + return State(θnp1,Ωnp1,x.t+p.Δt) +end + +function dsrun(x::State,p::Parameter,nstep::Int) + X = Vector{State}(undef,nstep) + for ii = 1:nstep + x = dsmap(x,p) + X[ii] = x + end + return X +end +``` + +## Simple pendulum +First we analyse the dynamics of the pendulum without anchor point movement, e.g. $a_0 = 0$ (in general) or $c = 0$ (for sinusoidal). We visualise the time series in terms of angle vs time. In this simplified case, we can easily visualise the entire state space, since it is two-dimensional, i.e. visualise angular velocity vs angle. + +```julia +using Plots +function plot_evolution_state_simple(X::Vector{State},p::Parameter) + t = getproperty.(X,:t) + θ = mod.((getproperty.(X,:θ).+π),2π).-π + Ω = getproperty.(X,:Ω) + p1 = plot( + t/τ,θ/π, + xlabel="t/τ",ylabel="θ/π", + seriestype=:scatter,marker=:circle,markersize=0.01, + legend=false + ) + p2 = plot( + θ/π,Ω*τ/π, + xlabel="θ/π",ylabel="Ωτ/π", + xlims=(-1.1,1.1),ylims=(-1.1,1.1), + seriestype=:scatter,marker=:circle,markersize=0.01, + aspect_ratio=1, + legend=false + ) + plot(p1,p2;layout=(1,2)) +end +``` + +### Periodic swing +We initialise the pendulum at a 90° angle with no initial velocity. This results in a periodic swinging motion which is very well captured by our numerical method. The state space is a circle whose centre is the stable equilibrium ($\theta=0, \Omega=0$). + +```julia +# Initial conditions +θ0 = π/2 +Ω0 = 0.0 +t0 = 0.0 +x = State(θ0,Ω0,t0) +# Parameter +τ = 1.0 +c = 0.0 +k = 0.0 +Δt = 0.01 +p = Parameter(τ,c,k,Δt) +# Simulate +nstep = 10000 +X = dsrun(x,p,nstep) +plot_evolution_state_simple(X,p) +``` + +### Unstable equilibrium +Next we initialise at the unstable equilibrium point ($\theta=\pi,\Omega=0$) and iterate in order to observe whether computational errors cause disturbances which lead to instability. The unstable state can be sustained for reasonable time using our numerical method and a sufficiently small time step. + +```julia +# Initial conditions +θ0 = π +Ω0 = 0.0 +t0 = 0.0 +x = State(θ0,Ω0,t0) +# Parameter +τ = 1.0 +c = 0.0 +k = 0.0 +Δt = 0.01 +p = Parameter(τ,c,k,Δt) +# Simulate +nstep = 10000 +X = dsrun(x,p,nstep) +plot_evolution_state_simple(X,p) +``` + +### Heteroclinics +Now we add a small disturbance to the unstable equilibrium. We observe the state of the system traveling along the heteroclinics between the stable and unstable equilibrium points. + +```julia +# Initial conditions +θ0 = π*(1-1e-6) +Ω0 = 0.0 +t0 = 0.0 +x = State(θ0,Ω0,t0) +# Parameter +τ = 1.0 +c = 0.0 +k = 0.0 +Δt = 0.01 +p = Parameter(τ,c,k,Δt) +# Simulate +nstep = 10000 +X = dsrun(x,p,nstep) +plot_evolution_state_simple(X,p) +``` + +## Pendulum with anchor point movement + +```julia +# Initial conditions +θ0 = π/2 +Ω0 = 0.0 +t0 = 0.0 +x = State(θ0,Ω0,t0) +# Parameter +τ = 1.0 +c = 0.75 +k = 1.0 +Δt = 0.01 +p = Parameter(τ,c,k,Δt) +# Simulate +nstep = 100000 +X = dsrun(x,p,nstep) +plot_evolution_state_simple(X,p) +```