### 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
# ╔═╡ 7ca39625-cba2-4637-b30c-5a1400216f12
begin
import Pkg
Pkg.activate("/home/hisacro/.julia/environments/v1.9/Project.toml")
using Plots, PlutoUI, IntervalRootFinding, IntervalArithmetic, StaticArrays, ImplicitPlots,ImplicitEquations, IntervalConstraintProgramming, FFTW, ModelingToolkit; plotly();
end
# ╔═╡ 7178185c-990e-11ee-3cd5-85a5c6b2f939
md"""
# Abstract
1. Added parameter α to account for rotors's leg length and re-doing the same analysis for fixed points as previous and in addition analyze $\frac{dr}{dt}=0$ cases
2. Introduced $(a-b\Delta \dot{\theta})$ into potential (Later this approach was found out to be redunact so a is set as 1, and b is set as 0)
# Goal
Finding points of the dynamical system.
# Methodology
1. Set all the governing ODEs simulatenously to 0.
2. Solve the equations, check for analytical solutions.
3. Divde the problem into multiple cases if present.
4. Categorize the fixed points of one case into trivial & non-trivial
"""
# ╔═╡ 627d1441-536d-45c0-99bc-a7e28f4774af
TableOfContents()
# ╔═╡ 5a68e070-3861-4e96-b010-d43b559a4f5f
md"""
##### adding packages
"""
# ╔═╡ 7a695d49-ada1-4e07-b529-535b14489464
# ╔═╡ 6c55169f-ef31-4315-b194-3a9ada720154
md"""
# Parameters
### Potential order, n = ($(@bind n Scrubbable(1:1:5,default=3)))
Length of rotor arms α = $(@bind α Slider(.1:0.05:2.; default=1.0, show_value=true))
Length between 2 rotors L = $(@bind l PlutoUI.NumberField(0.5:0.1:10., default=1.5; ))
Coupling Parameter K= $(@bind K PlutoUI.Slider(0.5:0.1:20., default=2.0; show_value=true))
Additive Parameter ξ= $(@bind ξ PlutoUI.NumberField(0.00:0.1:10., default=0.0; ))
#### Velocity parameters
Parameter A = $(@bind A1 PlutoUI.NumberField(0.5:0.1:20., default=1; ))
Parameter B = $(@bind B1 PlutoUI.NumberField(0.0:0.1:20., default=0.0; ))
#### 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; ))
"""
# ╔═╡ ae9771ac-47f0-45aa-ac49-5090f63fa1e9
md"""
# Potential
$V = \frac{1-e^{-kr^n}}{r}(a-b\Delta\dot{\theta})$
and
$
\frac{\partial V}{\partial r}=-\frac{e^{-kr^n}}{r^{2}}(e^{kr^n}-nKr^n-1)$
"""
# ╔═╡ f3dd1bdb-2d31-4bb9-8024-8c954859ef3b
md"""
# Potential Plots
"""
# ╔═╡ bcb38d85-9687-4c80-89a6-c89166f66e41
md"""
# Solutions of Trancendental Equations
"""
# ╔═╡ 79b0ea87-7091-4f0f-b9b5-d918fa6da2fb
md"""
r is now defined as
$r^2={\alpha^{2}+1+L^{2}-2L(\alpha cos\theta_{1}-cos\theta_{2})-2\alpha cos(\theta_{1}-\theta_{2})}$
"""
# ╔═╡ e2e1eb49-28a0-4979-a824-22a47790b0f1
begin
r(θ1,θ2)=sqrt(abs(α^2+1+l^2-2*l*(α*cos(θ1)-cos(θ2))-2*α*cos(θ1-θ2)));
end
# ╔═╡ c37c5c5d-0bc4-4045-afd7-affdb53adae0
begin
o1=range(-pi, stop=pi, length=500)
o2=range(-pi, stop=pi, length=500)
potential(θ1,θ2)=(1-exp(-K*(r(θ1,θ2)+ξ)^n))/(r(θ1,θ2)+ξ)
# potential(θ1,θ2)=(exp(-K*(r(θ1,θ2)+ξ)))/(r(θ1,θ2)+ξ)
end
# ╔═╡ 10b3ecb0-57e1-44d8-93c6-27e5701ec1ec
begin
Plots.plot(xlabel="θ1",ylabel="θ2",zlabel="potential")
# Plots.plot!(title="K=$(K), α=$(α), ξ=$(ξ), L=$(l), A=$(A1), B=$(B1)")
# Plots.plot!(title="order $(n), L=$(l), K=$(K), α=$(α)")
Plots.plot!(title="order $(n), L=$(l), K=$(K)")
# Plots.plot!(title="L=$(l), K=$(K)")
surface!(o1, o2, potential)
end
# ╔═╡ 48ac5e62-7d40-4f5a-8e95-546920ec7850
md"""
solutions of r are
$\begin{equation}
cos^{-1}\left(\frac{L^{2}+(\alpha^{2}-1)}{2L\alpha}\right),cos^{-1}\left(-\frac{L^{2}-(\alpha^{2}-1)}{2L}\right)
\end{equation}$
"""
# ╔═╡ 49d80e66-3efb-4a34-b225-f45bd5dc7e6d
begin
c=[acos( (l^2+α^2-1)/(2*l*α) ),acos((-l^2+α^2-1)/(2*l))]
@show c[1],c[2]
@show r(c[1],c[2])
end
# ╔═╡ dad92584-917a-4d51-aed2-1914f5952212
md"""
## Validity of solutions
"""
# ╔═╡ 91c0b0cf-7d54-4559-b857-2f3378261d0e
md""" we have relations
$αcos(θ_1-θ_2) = \frac{-L^2+α^2+1}{2}$
Subsituting Analytically found solutions into above, LHS = $(α*cos(c[1]-c[2])), RHS = $((-l^2+α^2+1)/2) and
$αcosθ_1-cosθ_2 = L$
LHS = $(α*cos(c[1])-cos(c[2]))
"""
# ╔═╡ 234fa59e-1a9c-4b98-b7e0-c9adf1a056ce
md"""
# Fixed points by linearization
"""
# ╔═╡ 93e3ffaf-28bc-49f2-a724-4efdfb546c30
md"""
ODEs governing equation of motion are
$
\begin{align*}
\dot{\theta}_{1} & =\frac{p_{1}}{\alpha^{2}}\\
\dot{p_{1}} & =-\frac{\partial V}{\partial r}\big(\frac{\partial r}{\partial\theta_{1}}(a-b\Delta\dot{\theta})+\dot{r}b)\\
\dot{\theta}_{2} & =p_{2}\\
\dot{p_{2}} & =-\frac{\partial V}{\partial r}\big(\frac{\partial r}{\partial\theta_{2}}(a-b\Delta\dot{\theta})-\dot{r}b)
\end{align*}$
and
$
\begin{align*}
\frac{\partial r}{\partial\theta_{1}} & =\alpha(Lsin\theta_{1}+sin(\theta_{1}-\theta_{2}))\frac{1}{r}\\
\frac{\partial r}{\partial\theta_{2}} & =(-Lsin\theta_{2}+\alpha sin(\theta_{2}-\theta_{1}))\frac{1}{r}\\
\dot{r} & =-\frac{1}{r}(-L\alpha sin\theta_{1}\dot{\theta_{1}}+Lsin\theta_{2}\dot{\theta_{2}})+\alpha sin(\theta_{1}-\theta_{2})(\dot{\theta_{1}}-\dot{\theta_{2}})
\end{align*}$
"""
# ╔═╡ f5e4f1f0-cfa7-4f21-adee-ec9efa5dc2db
md"""
## Case I: When trigonometric factors equal $0$
"""
# ╔═╡ 99cd4184-8ee3-499f-8e6c-435f6110fba5
md"""
### Numerically Approximated Roots of $\dot{p_1}$, $\dot{p_2}$
"""
# ╔═╡ 3b5dfb6c-c999-4983-8e68-732e4f6ad935
md"""
$\begin{align*}
\dot{p_{1}} = 0 & =\alpha(Lsin\theta_{1}+sin(\theta_{1}-\theta_{2}))\\
\dot{p_{2}} = 0 & = L sin\theta_{2}+\alpha sin(\theta_{2}-\theta_{1})
\end{align*}$
"""
# ╔═╡ 6f2a0a46-60e6-4b38-8e37-7554c287df87
function trans(z)
θ1,θ2 = z
SVector(sin(θ1-θ2)+l*sin(θ1), α*sin(θ2-θ1)-l*sin(θ2));
end
# ╔═╡ 3764bc0c-4c65-4956-b55e-f8745b5ba397
rts = IntervalRootFinding.roots(trans, IntervalBox(-pi..pi, 2));#,Newton, 1e-15)
# ╔═╡ 449ff82d-65e2-418c-a7b7-fcafab9b660e
md"""
Interval of Unique roots
"""
# ╔═╡ 87adb637-efa7-4d80-b393-e3f802492b58
[round.(interval(u),digits=10) for u ∈ rts if u.status == :unique]
# ╔═╡ 67bc0c4b-38f9-4b9c-b191-37bf70ad37d9
md"""
## Case II: when exponential factor go to $0$
"""
# ╔═╡ 7ee2f102-7ee8-435c-b61c-7712d3dd4337
md"""
This is just an abstract idea, it's yet to be developed. The potential taken here is the simpler version of Yukawa potential $V_0=\frac{e^-Kr}{r}$
"""
# ╔═╡ 40b932a8-f877-4236-a3aa-a5bfdb73a173
md"""
ODEs governing equation of motion were
$
\begin{align*}
\dot{\theta}_{1} & =\frac{p_{1}}{\alpha^{2}}\\
\dot{p_{1}} & =-\frac{\partial V_{0}}{\partial r}\big(\frac{\partial r}{\partial\theta_{1}}(a-b\Delta\dot{\theta})+\dot{r}b)\\
\dot{\theta}_{2} & =p_{2}\\
\dot{p_{2}} & =-\frac{\partial V_{0}}{\partial r}\big(\frac{\partial r}{\partial\theta_{2}}(a-b\Delta\dot{\theta})-\dot{r}b)
\end{align*}$
expanding out $\frac{\partial V_{0}}{\partial r}=-\frac{e^{-kr}}{r^{2}}(1+Kr)$,
$
\begin{align*}
\dot{\theta}_{1} & =\frac{p_{1}}{\alpha^{2}}\\
\dot{p_{1}} & =\frac{e^{-kr}}{r^{2}}(1+Kr)\big(\frac{\partial r}{\partial\theta_{1}}(a-b\Delta\dot{\theta})+\dot{r}b)\\
\dot{\theta}_{2} & =p_{2}\\
\dot{p_{2}} & =\frac{e^{-kr}}{r^{2}}(1+Kr)\big(\frac{\partial r}{\partial\theta_{2}}(a-b\Delta\dot{\theta})-\dot{r}b)
\end{align*}$
Taylor expand $e^{-Kr}$ and neglecting higher order terms,
$
\begin{align*}
\dot{\theta}_{1} & =\frac{p_{1}}{\alpha^{2}}\\
\dot{p_{1}} & =\frac{1}{r^{2}}(1-Kr)(1+Kr)\big(\frac{\partial r}{\partial\theta_{1}}(a-b\Delta\dot{\theta})+\dot{r}b)\\
\dot{\theta}_{2} & =p_{2}\\
\dot{p_{2}} & =\frac{1}{r^{2}}(1-Kr)(1+Kr)\big(\frac{\partial r}{\partial\theta_{2}}(a-b\Delta\dot{\theta})-\dot{r}b)
\end{align*}$
"""
# ╔═╡ e6e9894c-5637-4804-9954-b3c54b6168de
md"""
all time derivatives simulatenously go to zero, when $\dot{p_{1}}=\dot{p_{2}}=0$ and
$\begin{align}
0 & =1-Kr\\
r & =\frac{1}{K}\\
\end{align}$
"""
# ╔═╡ 02639bfc-aad5-47a4-809f-b381b3549d6b
md"""
### Solutions of $\frac{dr}{dt} = 0$
"""
# ╔═╡ df026d31-0bce-40e4-a7bd-e91d16626c18
md"""
$\begin{align*}
|\vec{r}|^{2} & =\alpha^{2}+1+L^{2}-2L(\alpha cos\theta_{1}-cos\theta_{2})-2\alpha cos(\theta_{1}-\theta_{2})\\
2r\dot{r} & =-2L(-\alpha sin\theta_{1}\dot{\theta_{1}}+sin\theta_{2}\dot{\theta_{2}})+2\alpha sin(\theta_{1}-\theta_{2})(\dot{\theta_{1}}-\dot{\theta_{2}})\\
0 & =-2L(-\alpha sin\theta_{1}\dot{\theta_{1}}+sin\theta_{2}\dot{\theta_{2}})+2\alpha sin(\theta_{1}-\theta_{2})(\dot{\theta_{1}}-\dot{\theta_{2}})
\end{align*}$
"""
# ╔═╡ 8f62179e-2d03-46bc-91a6-e7b5d69a1545
begin
r1(θ1,θ2)=sqrt((α^2+1+l^2-2*l*(α*cos(θ1)-cos(θ2))-2*α*cos(θ1-θ2)))
rg(θ1,θ2) = exp(K*r1(θ1,θ2)^n)-n*K*r1(θ1,θ2)^n-1
fg(θ1,θ2) = 2*l*(-α*sin(θ1)*ω1+sin(θ2)*ω2)+ 2*α*sin(θ1-θ2)*(ω1-ω2)
g1(θ1,θ2) = -l*(-sin(θ1))+sin(θ1-θ2)
g2(θ1,θ2) = -l*(sin(θ2))-α*sin(θ1-θ2)
Plots.plot(xlabel="θ1",ylabel="θ2")
# Plots.plot!(title=" L=$(l),α=$(α), K=$(K)")
# implicit_plot!(fg; xlim=(-pi,pi), ylim=(-pi,pi), lc=:purple, lw=2, ls=:solid, lab=:"dr/dt") #dr/dt
implicit_plot!(g1, xlim=(-pi,pi), ylim=(-pi,pi), ls=:dot, lw=2, lab=:"∂r/∂θ1") #dr/dθ_1 = ̇p_1
implicit_plot!(g2, xlim=(-pi,pi), ylim=(-pi,pi), lc=:red, lw=2, lab=:"∂r/∂θ2" ) #dr/dθ_2 = ̇p_2
implicit_plot!(rg, xlim=(-pi,pi), ylim=(-pi,pi), lc=:black, ls=:dashdot, lw=3, lab=:"∂v/∂r" )
end
# ╔═╡ eba7a32b-24e7-4e2b-9b3d-7023a06dc43f
begin
f(θ1,θ2) = α*(sin(θ1-θ2)+l*sin(θ1));
g(θ1,θ2) = α*sin(θ2-θ1)-l*sin(θ2);
Plots.plot(xlabel="θ1",ylabel="θ2")
# Plots.plot!(title=" L=$(l),α=$(α)")
# Plots.plot!(title=" L=$(l)")
implicit_plot!(g1, xlim=(-pi,pi), ylim=(-pi,pi), lw=:2, ls=:dot, lc=:blue, lab=:"dr/dθ1") #dr/dθ_1 = ̇p_1
implicit_plot!(g2, xlim=(-pi,pi), ylim=(-pi,pi), lw=:2, lc=:red, lab=:"dr/dθ2") #dr/dθ_2 = ̇p_2
# Plots.scatter!([c[1],-c[1]],[c[2],-c[2]],color="black", label=false)
# Plots.annotate!([(c[1]+0.25,c[2]-0.25,("S4",10))])
# Plots.annotate!([(-c[1]+0.25,-c[2]+0.25,("S4",10))])
Plots.scatter!([0,pi,pi,-pi,-pi],[0,pi,-pi,pi,-pi],color="black", label=false)
Plots.annotate!([(0.25,0,("S1",10))])
Plots.annotate!([(pi-0.25,pi+0,("S1",8))])
Plots.annotate!([(pi-0.25,-pi+0,("S1",8))])
Plots.annotate!([(-pi+0.25,pi+0,("S1",8))])
Plots.annotate!([(-pi+0.25,-pi+0,("S1",8))])
Plots.scatter!([pi,-pi],[0,0],color="black", label=false)
Plots.annotate!([(pi,0+0.25,("S2",10))])
Plots.annotate!([(-pi,0-0.25,("S2",10))])
Plots.scatter!([0,0],[pi,-pi],color="black", label=false)
Plots.annotate!([(0-0.3,pi-0.,("S3",10))])
Plots.annotate!([(0+0.3,-pi+0.0,("S3",10))])
Plots.plot!(xlim=(-pi-.1,pi+.1),ylim=(-pi-.1,pi+.1))
Plots.plot!(showaxis=:true,legend=false, grid=:false)
plot!(xtick=(-3.141592653589793:3.141592653589793:3.141592653589793,["-π",0,"π"]))
plot!(ytick=(-3.141592653589793:3.141592653589793:3.141592653589793,["-π",0,"π"]))
#plot!(ytick=pitick(theta2_start*pi, theta2_end*pi, 1))
# Plots.plot!(Eq(g, 0.), xlim=(-pi,pi), ylim=(-pi,pi),fc=:red)
end
# ╔═╡ 0aa5dbe6-d94a-45fc-af74-906a6e5688ad
# ╠═╡ disabled = true
#=╠═╡
begin
vars = @variables x y
#C1 = @constraint 1*(sin(x-y)+1.5*sin(x)) <= 0
#C2 = @constraint 1sin(y-x)-1.5*sin(y) <= 0
C1(x,y) = α*(sin(x-y)+l*sin(x)) < 0
C2(x,y) = α*sin(y-x)-l*sin(y) < 0
S1 = Separator(vars, C1)
S1 = Separator(vars, C2)
X = IntervalBox(-pi..pi, 2)
paving1 = pave(S1, X, 0.1)
paving2 = pave(S2, X, 0.1)
plot(paving1.boundary, lw=.01, lc=:blue, alpha=1, label="")
plot!(paving2.boundary, lw=.01, lc=:red, alpha=1, label="")
#plot!(paving2.inner)
end
╠═╡ =#
# ╔═╡ 1e8e50ea-d4f7-4421-97d5-f1c9659021da
# ╠═╡ disabled = true
# ╠═╡ skip_as_script = true
#=╠═╡
begin
function trans1(z)
θ1,θ2 = z
SVector(-2*l*(-α*sin(θ1)*ω1+sin(θ2)*ω2)+ 2*α*sin(θ1-θ2)*(ω1-ω2))
end
v = interval(0,1)
d = interval(1, 2)
box= IntervalBox(v,d)
rts1 = roots(trans1, box)
end
╠═╡ =#
# ╔═╡ df62a5c5-daf2-4ed5-b428-0a4648288693
# ╠═╡ disabled = true
#=╠═╡
@bind asdasd let
p = PlutoPlot(PlutoPlotly.Plot(PlutoPlotly.scatter(y = rand(10), name = "test", showlegend=true)))
add_plotly_listener!(p,"plotly_click", "
(e) => {
console.log(e)
let dt = e.points[0]
PLOT.value = [dt.x, dt.y]
PLOT.dispatchEvent(new CustomEvent('input'))
}
")
p
end
╠═╡ =#
# ╔═╡ 8189573b-f038-4b87-9385-bb62699ac16e
# ╠═╡ disabled = true
#=╠═╡
begin
data = 1;
p = plot([1,1.4,2,3],[1,.8,1,8])
on(p["click"]) do data
end
end
╠═╡ =#
# ╔═╡ ef2b261e-3129-45d9-b967-04815f597d77
md"""
# Force Plots
"""
# ╔═╡ c633d15f-9624-423c-bc1a-f53577b54615
begin
fx=range(0, stop=2, length=100)
force(r)=(exp(-K*r^n)*(exp(K*r^n)-K*n*(r^n)-1))/r^2
# force(r)=(K*r+1)*exp(-K*r)/r^2
fy=force.(fx)
Plots.plot(xlabel="r",ylabel="-∂V/∂r")
Plots.plot!(fx,fy,lab="force")
end
# ╔═╡ b4b85a7c-55e5-4cf2-9ac2-60e4f0b10727
md"""
velocity difference = $(round(ω1-ω2;sigdigits=4))
"""
# ╔═╡ Cell order:
# ╟─7178185c-990e-11ee-3cd5-85a5c6b2f939
# ╟─627d1441-536d-45c0-99bc-a7e28f4774af
# ╠═5a68e070-3861-4e96-b010-d43b559a4f5f
# ╠═7a695d49-ada1-4e07-b529-535b14489464
# ╠═7ca39625-cba2-4637-b30c-5a1400216f12
# ╟─6c55169f-ef31-4315-b194-3a9ada720154
# ╟─ae9771ac-47f0-45aa-ac49-5090f63fa1e9
# ╟─f3dd1bdb-2d31-4bb9-8024-8c954859ef3b
# ╟─c37c5c5d-0bc4-4045-afd7-affdb53adae0
# ╟─10b3ecb0-57e1-44d8-93c6-27e5701ec1ec
# ╟─bcb38d85-9687-4c80-89a6-c89166f66e41
# ╟─79b0ea87-7091-4f0f-b9b5-d918fa6da2fb
# ╟─e2e1eb49-28a0-4979-a824-22a47790b0f1
# ╟─48ac5e62-7d40-4f5a-8e95-546920ec7850
# ╟─49d80e66-3efb-4a34-b225-f45bd5dc7e6d
# ╟─dad92584-917a-4d51-aed2-1914f5952212
# ╟─91c0b0cf-7d54-4559-b857-2f3378261d0e
# ╟─234fa59e-1a9c-4b98-b7e0-c9adf1a056ce
# ╟─93e3ffaf-28bc-49f2-a724-4efdfb546c30
# ╟─f5e4f1f0-cfa7-4f21-adee-ec9efa5dc2db
# ╟─99cd4184-8ee3-499f-8e6c-435f6110fba5
# ╟─3b5dfb6c-c999-4983-8e68-732e4f6ad935
# ╟─eba7a32b-24e7-4e2b-9b3d-7023a06dc43f
# ╟─6f2a0a46-60e6-4b38-8e37-7554c287df87
# ╟─3764bc0c-4c65-4956-b55e-f8745b5ba397
# ╟─449ff82d-65e2-418c-a7b7-fcafab9b660e
# ╟─87adb637-efa7-4d80-b393-e3f802492b58
# ╟─67bc0c4b-38f9-4b9c-b191-37bf70ad37d9
# ╟─7ee2f102-7ee8-435c-b61c-7712d3dd4337
# ╟─40b932a8-f877-4236-a3aa-a5bfdb73a173
# ╟─e6e9894c-5637-4804-9954-b3c54b6168de
# ╟─02639bfc-aad5-47a4-809f-b381b3549d6b
# ╟─df026d31-0bce-40e4-a7bd-e91d16626c18
# ╟─8f62179e-2d03-46bc-91a6-e7b5d69a1545
# ╟─0aa5dbe6-d94a-45fc-af74-906a6e5688ad
# ╟─1e8e50ea-d4f7-4421-97d5-f1c9659021da
# ╠═df62a5c5-daf2-4ed5-b428-0a4648288693
# ╠═8189573b-f038-4b87-9385-bb62699ac16e
# ╟─ef2b261e-3129-45d9-b967-04815f597d77
# ╟─c633d15f-9624-423c-bc1a-f53577b54615
# ╟─b4b85a7c-55e5-4cf2-9ac2-60e4f0b10727