first implementation of pendulum
This commit is contained in:
Michael Krayer
parent
3bee6eb435
commit
1618df638f
|
|
@ -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)}$$
|
||||
<!--$$ c := \frac{c^{*}T^2}{l^{*}} $$-->
|
||||
$$ 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)
|
||||
```
|
||||
Loading…
Reference in New Issue