WARNING: Method definition default_nlsolve(Nothing, Base.Val{true}, Any, SciMLBase.NonlinearProblem{uType, isinplace, P, F, K, PT} where PT where K where F where P where isinplace where uType) in module OrdinaryDiffEqCore at /home/runner/.julia/packages/OrdinaryDiffEqCore/rRqPh/src/initialize_dae.jl:123 overwritten in module OrdinaryDiffEqNonlinearSolve at /home/runner/.julia/packages/OrdinaryDiffEqNonlinearSolve/0gbhD/src/initialize_dae.jl:1.
+ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
+WARNING: Method definition default_nlsolve(Nothing, Base.Val{true}, Any, SciMLBase.NonlinearProblem{uType, isinplace, P, F, K, PT} where PT where K where F where P where isinplace where uType) in module OrdinaryDiffEqCore at /home/runner/.julia/packages/OrdinaryDiffEqCore/rRqPh/src/initialize_dae.jl:123 overwritten in module OrdinaryDiffEqNonlinearSolve at /home/runner/.julia/packages/OrdinaryDiffEqNonlinearSolve/0gbhD/src/initialize_dae.jl:1.
+ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
+┌ Warning: Module OrdinaryDiffEqNonlinearSolve with build ID ffffffff-ffff-ffff-0000-008ddbad0360 is missing from the cache.
+│ This may mean OrdinaryDiffEqNonlinearSolve [127b3ac7-2247-4354-8eb6-78cf4e7c58e8] does not support precompilation but is imported by a module that does.
+└ @ Base loading.jl:1948
+┌ Warning: Module OrdinaryDiffEqNonlinearSolve with build ID ffffffff-ffff-ffff-0000-008ddbad0360 is missing from the cache.
+│ This may mean OrdinaryDiffEqNonlinearSolve [127b3ac7-2247-4354-8eb6-78cf4e7c58e8] does not support precompilation but is imported by a module that does.
+└ @ Base loading.jl:1948
+┌ Warning: Module OrdinaryDiffEqSDIRK with build ID ffffffff-ffff-ffff-0000-008f70ed8c78 is missing from the cache.
+│ This may mean OrdinaryDiffEqSDIRK [2d112036-d095-4a1e-ab9a-08536f3ecdbf] does not support precompilation but is imported by a module that does.
+└ @ Base loading.jl:1948
+┌ Warning: Module OrdinaryDiffEqBDF with build ID ffffffff-ffff-ffff-0000-008ff8fdd020 is missing from the cache.
+│ This may mean OrdinaryDiffEqBDF [6ad6398a-0878-4a85-9266-38940aa047c8] does not support precompilation but is imported by a module that does.
+└ @ Base loading.jl:1948
+┌ Warning: Module OrdinaryDiffEqNonlinearSolve with build ID ffffffff-ffff-ffff-0000-008ddbad0360 is missing from the cache.
+│ This may mean OrdinaryDiffEqNonlinearSolve [127b3ac7-2247-4354-8eb6-78cf4e7c58e8] does not support precompilation but is imported by a module that does.
+└ @ Base loading.jl:1948
+┌ Warning: Module OrdinaryDiffEqNonlinearSolve with build ID ffffffff-ffff-ffff-0000-008ddbad0360 is missing from the cache.
+│ This may mean OrdinaryDiffEqNonlinearSolve [127b3ac7-2247-4354-8eb6-78cf4e7c58e8] does not support precompilation but is imported by a module that does.
+└ @ Base loading.jl:1948
+┌ Warning: Module OrdinaryDiffEqNonlinearSolve with build ID ffffffff-ffff-ffff-0000-008ddbad0360 is missing from the cache.
+│ This may mean OrdinaryDiffEqNonlinearSolve [127b3ac7-2247-4354-8eb6-78cf4e7c58e8] does not support precompilation but is imported by a module that does.
+└ @ Base loading.jl:1948
+┌ Warning: Module OrdinaryDiffEq with build ID ffffffff-ffff-ffff-0000-008cb1f2f77f is missing from the cache.
+│ This may mean OrdinaryDiffEq [1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] does not support precompilation but is imported by a module that does.
+└ @ Base loading.jl:1948Let us plot some extremals, solutions of this flow. The initial condition $x_0$ is fixed while we compute some extremals for different values of initial covector $p_0$. We compute some specific initial covectors for a nice plot.
Here, the shooting equation given by
\[ S({p₀}) = \pi(z(t_f,x_0,{p₀})) - x_f = 0,\]
with $\pi(x, p) = x$, has two solutions: $p₀ = -0.9851$ and $p₀ = 0.5126$.
π((x, p)) = x
+plot(plt_x, plt_p, plt_u, plt_phase; layout=(2, 2), size=(800, 600))
Here, the shooting equation given by
\[ S({p₀}) = \pi(z(t_f,x_0,{p₀})) - x_f = 0,\]
with $\pi(x, p) = x$, has two solutions: $p₀ = -0.9851$ and $p₀ = 0.5126$.
π((x, p)) = x
# Shooting function
S(p0) = (π ∘ ocp_flow)(t0, x0, p0, tf) - xf
@@ -169,7 +196,7 @@
plot!(plt2_u; xlabel="t", ylabel="u(t,p₀)", legend=false, ylims=(-6, 5))
plot!(plt2_phase; xlabel="x(t,p₀)", ylabel="p(t,p₀)", legend=false, xlims=(0, 2.5), ylims=(-1, 5))
-plot(plt2_x, plt2_p, plt2_u, plt2_phase; layout=(2, 2), size=(800, 600))
Now, we can compute the conjugate points along the two extremals. We have to compute the flow $\delta z(t, p₀)$ of the Jacobi equation with the initial condition $\delta z(0) = (0, 1)$. This is given solving
\[ \delta z(t, p₀) = \dfrac{\partial}{\partial p₀}z(t, p₀).\]
Note that to compute the conjugate points, we only need the first component:
\[ \delta z(t, p₀)_1.\]
using ForwardDiff
+plot(plt2_x, plt2_p, plt2_u, plt2_phase; layout=(2, 2), size=(800, 600))
Now, we can compute the conjugate points along the two extremals. We have to compute the flow $\delta z(t, p₀)$ of the Jacobi equation with the initial condition $\delta z(0) = (0, 1)$. This is given solving
\[ \delta z(t, p₀) = \dfrac{\partial}{\partial p₀}z(t, p₀).\]
Note that to compute the conjugate points, we only need the first component:
\[ \delta z(t, p₀)_1.\]
using ForwardDiff
function jacobi_flow(t, p0)
x(t, p0) = (π ∘ ocp_flow)(t0, x0, p0, t)
@@ -198,11 +225,11 @@
#
plt_conj = plot(plt_conj1, plt_conj2;
- layout=(1, 2), size=(800, 300), leftmargin=25px, bottommargin=15px)
We compute the first conjugate point along the first extremal and add it to the plot.
tau0 = Roots.find_zero(tau -> jacobi_flow(tau, sol1_p0), (0.4, 0.6))
+ layout=(1, 2), size=(800, 300), leftmargin=25px, bottommargin=15px)
We compute the first conjugate point along the first extremal and add it to the plot.
tau0 = Roots.find_zero(tau -> jacobi_flow(tau, sol1_p0), (0.4, 0.6))
println("For p0 = ", sol1_p0, " tau_0 = ", tau0)
-plot!(plt_conj[1], [tau0], [jacobi_flow(tau0, sol1_p0)]; seriestype=:scatter)
To conclude on this example, we compute the conjugate locus by using a path following algorithm. Define $F(\tau,p₀) = \delta x(\tau,p₀)$ and suppose that the partial derivative $\partial_\tau F(\tau,p₀)$ is invertible, then, by the implicit function theorem the conjugate time is a function of $p₀$. So, since here $p₀\in\R$, we can compute them by solving the initial value problem for $p₀ \in [\alpha, \beta]$:
\[ \dot{\tau}(p₀) = -\dfrac{\partial F}{\partial \tau}(\tau(p₀),p₀)^{-1}\,
+plot!(plt_conj[1], [tau0], [jacobi_flow(tau0, sol1_p0)]; seriestype=:scatter)
To conclude on this example, we compute the conjugate locus by using a path following algorithm. Define $F(\tau,p₀) = \delta x(\tau,p₀)$ and suppose that the partial derivative $\partial_\tau F(\tau,p₀)$ is invertible, then, by the implicit function theorem the conjugate time is a function of $p₀$. So, since here $p₀\in\R$, we can compute them by solving the initial value problem for $p₀ \in [\alpha, \beta]$:
\[ \dot{\tau}(p₀) = -\dfrac{\partial F}{\partial \tau}(\tau(p₀),p₀)^{-1}\, \dfrac{\partial F}{\partial p₀}(\tau(p₀),p₀), \quad \tau(\alpha) = \tau_0.\]
For the numerical experiment, we set $\alpha = -0.9995$, $\beta = -0.5$.
function conjugate_times_rhs_path(tau, p0)
dF = ForwardDiff.gradient(y -> jacobi_flow(y...), [tau, p0])
@@ -236,4 +263,4 @@
#
plot(plt_x, plt_conj_times;
- layout=(1,2), legend=false, size=(800,300), leftmargin=25px, bottommargin=15px)
Settings
This document was generated with Documenter.jl version 1.7.0 on Friday 6 September 2024. Using Julia version 1.10.5.