### A Pluto.jl notebook ###
# v0.19.35
using Markdown
using InteractiveUtils
# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
quote
local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
local el = $(esc(element))
global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
el
end
end
# ╔═╡ 4981e181-25d3-41cb-9803-fdf1e41c6b92
begin
import Pkg
Pkg.activate("/home/hisacro/.julia/environments/v1.9/Project.toml")
using Markdown
using InteractiveUtils
using ImplicitPlots, PlutoPlotly, HypertextLiteral
using ImplicitPlots.Contour
using PlutoUI
using ForwardDiff, LinearAlgebra
using DelimitedFiles, DynamicalSystems, OrdinaryDiffEq, ProgressMeter, FFTW
import Plots.scatter! as sscater
import Plots.plot! as pplot
import Plots.plot as fplot
end
# ╔═╡ 8d689808-d40d-4267-87d7-b6284520fed5
function plot_implicit(f::ImplicitFunction{2}; resolution = 400, lab = "", kwargs...)
# We keep this names but plotly does not use the name xlims or ylims for specifing axes limits. It uses the xaxis_range and yaxis_range ones
xlims = get(kwargs, :xlims, (-pi, pi))
ylims = get(kwargs, :ylims, xlims)
l = Layout(; template = "none", kwargs...)
l.xaxis_range = xlims
# We default to equal axis ratio and we set ylims
l.yaxis = attr(;
range = ylims,
scaleanchor = "x",
scaleratio = 1,
)
rx = range(xlims[1]; stop = xlims[2], length = resolution)
ry = range(ylims[1]; stop = ylims[2], length = resolution)
z = [f(x, y) for x in rx, y in ry]
lvl = Contour.contour(collect(rx), collect(ry), z, 0.0)
lines = Contour.lines(lvl)
clr = "dodgerblue"
data = map(enumerate(lines)) do (k, line)
xs, ys = Contour.coordinates(line)
scatter(;
line = attr(;width = 1, color = clr),
showlegend = k == 1,
name = k > 1 ? "" : lab,
x = xs,
y = ys,
)
end
plot(data,l)
end
# ╔═╡ c3325f5d-0ba5-4c27-9104-d23d357f5626
TableOfContents()
# ╔═╡ 0a998633-7070-46a3-a792-e8f8fd1a2715
md"""
# Parameters
### Potential order, n = ($(@bind n Scrubbable(1:1:5,default=3)))
### Exponential factor, K = $(@bind K PlutoUI.NumberField(1.5:0.1:20.0,default=2.0))
### Length between 2 rotors L = $(@bind l PlutoUI.NumberField(0.5:0.001:10., default=1.5; ))
### Additive parameter ξ = $(@bind ξ PlutoUI.NumberField(0.0:0.01:10., default=0.0; ))
#### Velocity parameters
Parameter A = $(@bind A1 PlutoUI.Slider(-1.0:0.1:20., default=1.0; show_value=true))
Parameter B = $(@bind B1 PlutoUI.Slider(0.0:0.1:20., default=0.0; show_value=true))
#### Velocities
$\dot\theta_1$ = $(@bind ω1 PlutoUI.NumberField(-10:0.05:10., default=2.05; ))
$\dot\theta_2$ = $(@bind ω2 PlutoUI.NumberField(-10.0:0.05:10., default=2.00; ))
#### α - rotor arm parameter = $(@bind α PlutoUI.Slider(0.2:0.01:1.0,default=1.0,show_value=true))
### For Evolution
Total Energy= $(@bind energy PlutoUI.NumberField(0.0:0.1:10., default=0.5; ))
Potential Energy= $(@bind p_energy PlutoUI.NumberField(0.0:0.001:2., default=0.5;))
p2 = $(@bind p2 PlutoUI.NumberField(-10.0:0.1:10., default=0.0; ))
"""
# ╔═╡ 50961da2-1a6f-4df7-901b-133cc5fcaf7c
md"""
# Goal
Finding valid initial conditions for the dynamical system.
# Methodology
## #1
1. Fix the energy $\mathcal{H}$ and then potential energy $V$ of the system such that $\mathcal{H}-V\geq0$.
2. Find $(\theta_{1},\theta_{2})$ corresponding to above conditions from $r(\theta_{1},\theta_{2})$ constrain.
3. Plot the contour of satisfing $(\theta_{1},\theta_{2})$, choose a single set using PlutoPlotly's point & click function.
4. Find $p_{1}$ by fixing $p_{2}$ ($(p2)),
$
p_{1}^{2}=2\alpha^{2}(\mathcal{H}-V-\frac{p_{2}^{2}}{2}).$
4. Iterate over $V$ and make a list of initial condtions, $(\theta_{1}, p_{1}, \theta_{2}, p_{2})$ for corresponding energy.
## #2
It has been mentioned in $\textcolor{red}{rotor\_transcend\_fixed.jl}$, the condition of $\frac{\partial V}{\partial r}=0$ leads to fixed points. This leads to an exotic system behaviour where there are infinite fixed points, due to limited scope
it's not explored further. Here for the sake of completion, such initial conditions are choosen.
1. Plot the contour $\frac{\partial V}{\partial r}=0$.
2. Choose any one of points that satisfy above (Point & click funtion with PlutoPlotly)
"""
# ╔═╡ a8a2555a-9d4d-4cd8-975a-35d42089fcfa
k=K;
# ╔═╡ 3a56e612-3306-4e3e-9020-9c7479b7acb4
md"""
# #1 Rotor plot
"""
# ╔═╡ 62baca48-1f41-4f1d-9fcd-3f63596bb65f
function polar2cart(Points)
y1 = α*sin.(Points[1])
x1 = α*cos.(Points[1])
y2 = sin.(Points[3])
x2 = cos.(Points[3])
((x1, y1), (l .+ x2, y2))
end
# ╔═╡ fcf5fecd-aaab-4d6b-8842-fafa185d97ab
md"""
## Potential Counters (Point & Click)
"""
# ╔═╡ ac759b72-88ae-49ce-9fbb-4bf9f439029c
let
f(x,y) = x^2+y^2-3
r(x,y)= sqrt((α^2+1+l^2-2*l*(α*cos(x)-cos(y))-2*α*cos(x-y)));
hv(x,y)=-p_energy+((1-exp(-K*(r(x,y)+ξ)^n))/(r(x,y)+ξ))
p = plot_implicit(ImplicitFunction(hv); xlims = (0, 2*pi), title="Potential Energy = $(p_energy)", lab = "K=$(k) L=$(l) α=$(α)", xaxis_title="θ_1", yaxis_title="θ_2", lw=2)
add_plotly_listener!(p, "plotly_click", "
(e) => {
let dt = e.points[0]
PLOT.value = [dt.x, dt.y]
PLOT.dispatchEvent(new CustomEvent('input'))
}
")
@bind point p
end
# ╔═╡ 770169d4-c7df-44f8-a8e6-df08de3497b7
point
# ╔═╡ 63e2774c-0ef8-420b-a56f-0fd33fde7579
begin
r(θ1,θ2)=sqrt(abs(α^2+1+l^2-2*l*(α*cos(θ1)-cos(θ2))-2*α*cos(θ1-θ2)));
end
# ╔═╡ 04769387-fe2d-4d2e-bd10-adee37f620df
md"""
#### Distance between 2 rotors at this point $(r(point[1],point[2]))
"""
# ╔═╡ c74db737-f073-4ddd-ad90-5a5a0b098dd5
md"""
# Initial Conditions
"""
# ╔═╡ 3028b6c4-9f32-4a80-99cb-481b989b272e
p1=sqrt(2*α^2*(energy-p_energy-((p2^2)/2)))
# ╔═╡ 0c37e832-34ff-4d09-bcd7-23c571ad4458
u0=[point[1],p1,point[2],p2]
# ╔═╡ b8ac45f0-7fc2-45b2-8cb1-97018bf8bd1e
rotor1_cord, rotor2_cord = polar2cart(u0)
# ╔═╡ 04ca4082-0300-4d35-bb7b-15e004e358ec
begin
rotor1_x, rotor1_y = rotor1_cord[1], rotor1_cord[2]
rotor2_x, rotor2_y = rotor2_cord[1], rotor2_cord[2]
end
# ╔═╡ 4fdc54bf-a090-4e22-ba96-0cc2d2dda4b0
begin
fplot(
([0,rotor1_x],[0,rotor1_y]),
color = [:black,:red],
xlim = (-5,5),
ylim = (-2,2)
)
sscater(
([0,rotor1_x],[0,rotor1_y]),
color = [:black,:red],
xlim = (-5,5),
ylim = (-2,2)
)
pplot(
([l,rotor2_x],[0,rotor2_y]),
color = [:black,:blue],
xlim = (-5,5),
ylim = (-2,2),
legend = false,
framestyle=:box
)
sscater(
([l,rotor2_x],[0,rotor2_y]),
color = [:black,:blue],
xlim = (-5,5),
ylim = (-2,2),
title ="Combined Rotor system"
)
end
# ╔═╡ 3282bbc8-bd09-41ef-90d1-837de5808915
@show u0
# ╔═╡ 25a82a6a-6103-4d0c-94da-f86281f89fb4
md"""
## Lyapunov exponent starting at these points
"""
# ╔═╡ 045a4ae6-f2e7-4b5b-8e41-9fc9663c28ff
md"""
## Time evolution parameters
"""
# ╔═╡ 0f0f74e3-1e7a-445e-b9db-57d54c3e7748
md"""
To monitor the changes in Lyapunov exponent over the period of evolution.
"""
# ╔═╡ 404e4156-4aad-4bb3-b38f-ec0e8837bf0d
function rotor_2(u, params, t)
θ1,p1,θ2,p2 = u
K,L,α,ξ,A,B = params
r=sqrt(abs(α^2+1+L^2-2*L*(α*cos(θ1)-cos(θ2))-2*α*cos(θ1-θ2)));
r=r+ξ
# v0 = (1-exp(-K*r^n))/r
pv0_r = -exp(-K*r^n)*(-1-n*K*r^n+exp(K*r^n))/(r^2)
pr_θ1 = (α*(sin(θ1-θ2)+L*sin(θ1)))/r
pr_θ2 = (α*sin(θ2-θ1) -L*sin(θ2)) /r
du1 = (p1)/(α^2)
du2 = -pv0_r*pr_θ1
du3 = (p2)
du4 = -pv0_r*pr_θ2
return SVector{4}(du1, du2, du3, du4)
end
# ╔═╡ 17992843-fc80-4f4a-9e8f-cf6689ea4e0a
params = [K,l,α,ξ,A1,B1]
# ╔═╡ 59cee79b-d7b7-414a-a798-cae8b9502f90
md"""
Time of Evolution $(@bind time_end PlutoUI.Slider(10e2:10e2:10e4,default=10e2,show_value=true))
Sampling time $(@bind sampling_time PlutoUI.NumberField(0.005:0.1:10., default=0.005; ))
Transient time $(@bind transient_time PlutoUI.Slider(0:10e2:10e4,default=0,show_value=true))
### inter = $(@bind inter PlutoUI.Slider(0:50:100000,default=100;show_value=true))
"""
# ╔═╡ 1a6be30e-8536-4fec-9649-7eeca2fb4acb
begin
algo = Vern7();
diffeq = (alg = algo, reltol = 1.0e-7, abstol = 1.0e-7, atol = 1.0e-7, rtol = 1.0e-7, )
times = (0, time_end)
#
end
# ╔═╡ 8064244c-9e1b-436e-b272-a1cfae3c6289
md"""
## rotor system
"""
# ╔═╡ 2fa79d10-3b4a-4526-9a92-bd96ad991372
rotor_system_1 = CoupledODEs(rotor_2, u0, params; diffeq, t0 = 0.0)
# ╔═╡ 3b1fbc92-1a8b-488b-9e2a-fc6a5af6026e
a=lyapunov(rotor_system_1,time_end; u0=u0, Ttr=transient_time, Δt=sampling_time)
# ╔═╡ 2b6167bf-667f-4de0-af65-0f1dc5c32d8e
md"""
Largest Lyapunov = $(a)
"""
# ╔═╡ 8b2c29c4-87c7-4b90-83e6-ddf13fd87edf
Points_1, t = trajectory(rotor_system_1, time_end; Ttr = transient_time, Δt =sampling_time)
# ╔═╡ 76929141-3476-4e96-9a74-b4685206b5f5
end_ind = lastindex(t)
# ╔═╡ 300976bb-842f-438f-92a6-110a5b0c99f6
ttt=plot(cos.(Points_1[:,1][1:inter:end_ind]), cos.(Points_1[:,3][1:inter:end_ind]), xlabel="cos(θ1)", ylabel="cos(θ2)", title="$(round.(u0,digits=5))", marker=1)
# ╔═╡ b2842351-0073-46ed-a3e4-4acc2348c6e7
md"""
# #2 Case of $\frac{\partial V}{\partial r}=0$ (Point & Click)
"""
# ╔═╡ 7ba12727-8238-4ef8-b491-93eafd3361bd
let
f(x,y) = x^2+y^2-3
rg(x,y) = exp(K*r(x,y)^n)-n*K*r(x,y)^n-1
p1 = plot_implicit(ImplicitFunction(rg); xlims = (-pi,pi), title=" case II", lab = "K=$(k) L=$(l) α=$(α)",xlabel="θ1",ylabel="θ2")
add_plotly_listener!(p1, "plotly_click", "
(e) => {
let dt = e.points[0]
PLOT.value = [dt.x, dt.y]
PLOT.dispatchEvent(new CustomEvent('input'))
}
")
@bind point_vr p1
end
# ╔═╡ 8ceae86f-3222-4014-a636-71b415508080
begin
u01 = [point_vr[1],0.0,point_vr[2],0.0]
@show u01
end
# ╔═╡ 38dc0db2-f159-41df-8257-da9043c1b137
# ╠═╡ disabled = true
#=╠═╡
J0_c1 = ForwardDiff.jacobian(f1_brute, u01)
╠═╡ =#
# ╔═╡ 8577e825-35c2-4fc7-9196-ed704c249e91
#=╠═╡
eigvals(J0_c1)
╠═╡ =#
# ╔═╡ Cell order:
# ╠═50961da2-1a6f-4df7-901b-133cc5fcaf7c
# ╟─4981e181-25d3-41cb-9803-fdf1e41c6b92
# ╟─8d689808-d40d-4267-87d7-b6284520fed5
# ╠═c3325f5d-0ba5-4c27-9104-d23d357f5626
# ╟─0a998633-7070-46a3-a792-e8f8fd1a2715
# ╟─a8a2555a-9d4d-4cd8-975a-35d42089fcfa
# ╟─3a56e612-3306-4e3e-9020-9c7479b7acb4
# ╟─62baca48-1f41-4f1d-9fcd-3f63596bb65f
# ╟─b8ac45f0-7fc2-45b2-8cb1-97018bf8bd1e
# ╟─04ca4082-0300-4d35-bb7b-15e004e358ec
# ╟─4fdc54bf-a090-4e22-ba96-0cc2d2dda4b0
# ╟─3282bbc8-bd09-41ef-90d1-837de5808915
# ╟─fcf5fecd-aaab-4d6b-8842-fafa185d97ab
# ╟─ac759b72-88ae-49ce-9fbb-4bf9f439029c
# ╟─04769387-fe2d-4d2e-bd10-adee37f620df
# ╟─770169d4-c7df-44f8-a8e6-df08de3497b7
# ╟─2b6167bf-667f-4de0-af65-0f1dc5c32d8e
# ╟─63e2774c-0ef8-420b-a56f-0fd33fde7579
# ╟─c74db737-f073-4ddd-ad90-5a5a0b098dd5
# ╠═3028b6c4-9f32-4a80-99cb-481b989b272e
# ╠═0c37e832-34ff-4d09-bcd7-23c571ad4458
# ╟─25a82a6a-6103-4d0c-94da-f86281f89fb4
# ╟─3b1fbc92-1a8b-488b-9e2a-fc6a5af6026e
# ╟─045a4ae6-f2e7-4b5b-8e41-9fc9663c28ff
# ╟─0f0f74e3-1e7a-445e-b9db-57d54c3e7748
# ╟─404e4156-4aad-4bb3-b38f-ec0e8837bf0d
# ╟─17992843-fc80-4f4a-9e8f-cf6689ea4e0a
# ╟─1a6be30e-8536-4fec-9649-7eeca2fb4acb
# ╟─76929141-3476-4e96-9a74-b4685206b5f5
# ╟─59cee79b-d7b7-414a-a798-cae8b9502f90
# ╟─8064244c-9e1b-436e-b272-a1cfae3c6289
# ╟─2fa79d10-3b4a-4526-9a92-bd96ad991372
# ╟─8b2c29c4-87c7-4b90-83e6-ddf13fd87edf
# ╠═300976bb-842f-438f-92a6-110a5b0c99f6
# ╟─b2842351-0073-46ed-a3e4-4acc2348c6e7
# ╠═7ba12727-8238-4ef8-b491-93eafd3361bd
# ╟─8ceae86f-3222-4014-a636-71b415508080
# ╟─38dc0db2-f159-41df-8257-da9043c1b137
# ╠═8577e825-35c2-4fc7-9196-ed704c249e91