@@ -7,9 +7,8 @@ const t0 = 0
77const tf = 1
88const x0 = - 1
99const xf = 0
10- const α = 1.5
10+ const α = 1.5
1111ocp = @def begin
12-
1312 t ∈ [t0, tf], time
1413 x ∈ R, state
1514 u ∈ R, control
@@ -19,8 +18,7 @@ ocp = @def begin
1918
2019 ẋ (t) == - x (t) + α * x (t)^ 2 + u (t)
2120
22- ∫ ( 0.5 u (t)^ 2 ) → min
23-
21+ ∫ (0.5 u (t)^ 2 ) → min
2422end
2523nothing # hide
2624
@@ -29,7 +27,7 @@ u(x, p) = p
2927
3028π ((x, p)) = x
3129
32- S (p0) = π ( φ (t0, x0, p0, tf) ) - xf # shooting function
30+ S (p0) = π (φ (t0, x0, p0, tf)) - xf # shooting function
3331
3432ξ = [0.1 ] # initial guess
3533
@@ -40,19 +38,20 @@ function fsolve(f, j, x; kwargs...)
4038 catch e
4139 println (" Erreur using MINPACK" )
4240 println (e)
43- println (" hybrj not supported. Replaced by hybrd even if it is not visible on the doc." )
41+ println (
42+ " hybrj not supported. Replaced by hybrd even if it is not visible on the doc."
43+ )
4444 MINPACK. fsolve (f, x; kwargs... )
4545 end
4646end
4747
4848using DifferentiationInterface
49- import ForwardDiff
49+ using ForwardDiff : ForwardDiff
5050backend = AutoForwardDiff ()
5151
52- nle! = ( s, ξ) -> s[1 ] = S (ξ[1 ]) # auxiliary function
52+ nle! = (s, ξ) -> s[1 ] = S (ξ[1 ]) # auxiliary function
5353jnle! = (js, ξ) -> jacobian! (nle!, similar (ξ), js, backend, ξ) #
5454
55-
5655indirect_sol = fsolve (nle!, jnle!, ξ; show_trace= true ) # resolution of S(p0) = 0
5756p0_sol = indirect_sol. x[1 ] # costate solution
5857println (" \n costate: p0 = " , p0_sol)
@@ -61,71 +60,85 @@ println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
6160sol = φ ((t0, tf), x0, p0_sol)
6261plot (sol)
6362
64- using Plots. PlotMeasures
65- function Plots. plot (S:: Function , p0:: Float64 ; Np0= 20 , kwargs... )
66-
63+ using Plots. PlotMeasures
64+ function Plots. plot (S:: Function , p0:: Float64 ; Np0= 20 , kwargs... )
65+
6766 # times for wavefronts
68- times = range (t0, tf, length= 3 )
67+ times = range (t0, tf; length= 3 )
6968
7069 # times for trajectories
71- tspan = range (t0, tf, length= 100 )
70+ tspan = range (t0, tf; length= 100 )
7271
7372 # interval of initial covector
74- p0_min = - 0.5
75- p0_max = 2
73+ p0_min = - 0.5
74+ p0_max = 2
7675
7776 # covector solution
78- p0_sol = p0
79-
77+ p0_sol = p0
78+
8079 # plot of the flow in phase space
81- plt_flow = plot ()
82- p0s = range (p0_min, p0_max, length= Np0)
83- for i ∈ eachindex (p0s)
80+ plt_flow = plot ()
81+ p0s = range (p0_min, p0_max; length= Np0)
82+ for i in eachindex (p0s)
8483 sol = φ ((t0, tf), x0, p0s[i])
8584 x = state (sol).(tspan)
8685 p = costate (sol).(tspan)
87- label = i== 1 ? " extremals" : false
88- plot! (plt_flow, x, p, color= :blue , label= label)
89- end
90-
86+ label = i== 1 ? " extremals" : false
87+ plot! (plt_flow, x, p; color= :blue , label= label)
88+ end
89+
9190 # plot of wavefronts in phase space
92- p0s = range (p0_min, p0_max, length= 200 )
93- xs = zeros (length (p0s), length (times))
94- ps = zeros (length (p0s), length (times))
95- for i ∈ eachindex (p0s)
96- sol = φ ((t0, tf), x0, p0s[i], saveat= times)
97- xs[i, :] .= state (sol).(times)
98- ps[i, :] .= costate (sol).(times)
99- end
100- for j ∈ eachindex (times)
101- label = j== 1 ? " flow at times" : false
102- plot! (plt_flow, xs[:, j], ps[:, j], color= :green , linewidth= 2 , label= label)
103- end
104-
91+ p0s = range (p0_min, p0_max; length= 200 )
92+ xs = zeros (length (p0s), length (times))
93+ ps = zeros (length (p0s), length (times))
94+ for i in eachindex (p0s)
95+ sol = φ ((t0, tf), x0, p0s[i]; saveat= times)
96+ xs[i, :] .= state (sol).(times)
97+ ps[i, :] .= costate (sol).(times)
98+ end
99+ for j in eachindex (times)
100+ label = j== 1 ? " flow at times" : false
101+ plot! (plt_flow, xs[:, j], ps[:, j]; color= :green , linewidth= 2 , label= label)
102+ end
103+
105104 #
106- plot! (plt_flow, xlims= (- 1.1 , 1 ), ylims= (p0_min, p0_max))
107- plot! (plt_flow, [0 , 0 ], [p0_min, p0_max], color= :black , xlabel= " x" , ylabel= " p" , label= " x=xf" )
108-
105+ plot! (plt_flow; xlims= (- 1.1 , 1 ), ylims= (p0_min, p0_max))
106+ plot! (
107+ plt_flow,
108+ [0 , 0 ],
109+ [p0_min, p0_max];
110+ color= :black ,
111+ xlabel= " x" ,
112+ ylabel= " p" ,
113+ label= " x=xf" ,
114+ )
115+
109116 # solution
110117 sol = φ ((t0, tf), x0, p0_sol)
111118 x = state (sol).(tspan)
112119 p = costate (sol).(tspan)
113- plot! (plt_flow, x, p, color= :red , linewidth= 2 , label= " extremal solution" )
114- plot! (plt_flow, [x[end ]], [p[end ]], seriestype= :scatter , color= :green , label= false )
115-
120+ plot! (plt_flow, x, p; color= :red , linewidth= 2 , label= " extremal solution" )
121+ plot! (plt_flow, [x[end ]], [p[end ]]; seriestype= :scatter , color= :green , label= false )
122+
116123 # plot of the shooting function
117- p0s = range (p0_min, p0_max, length= 200 )
118- plt_shoot = plot (xlims= (p0_min, p0_max), ylims= (- 2 , 4 ), xlabel= " p₀" , ylabel= " y" )
119- plot! (plt_shoot, p0s, S, linewidth= 2 , label= " S(p₀)" , color= :green )
120- plot! (plt_shoot, [p0_min, p0_max], [0 , 0 ], color= :black , label= " y=0" )
121- plot! (plt_shoot, [p0_sol, p0_sol], [- 2 , 0 ], color= :black , label= " p₀ solution" , linestyle= :dash )
122- plot! (plt_shoot, [p0_sol], [0 ], seriestype= :scatter , color= :green , label= false )
123-
124+ p0s = range (p0_min, p0_max; length= 200 )
125+ plt_shoot = plot (; xlims= (p0_min, p0_max), ylims= (- 2 , 4 ), xlabel= " p₀" , ylabel= " y" )
126+ plot! (plt_shoot, p0s, S; linewidth= 2 , label= " S(p₀)" , color= :green )
127+ plot! (plt_shoot, [p0_min, p0_max], [0 , 0 ]; color= :black , label= " y=0" )
128+ plot! (
129+ plt_shoot,
130+ [p0_sol, p0_sol],
131+ [- 2 , 0 ];
132+ color= :black ,
133+ label= " p₀ solution" ,
134+ linestyle= :dash ,
135+ )
136+ plot! (plt_shoot, [p0_sol], [0 ]; seriestype= :scatter , color= :green , label= false )
137+
124138 # final plot
125- plot (plt_flow, plt_shoot; layout= (1 ,2 ), leftmargin= 15 px, bottommargin= 15 px, kwargs... )
126-
139+ plot (plt_flow, plt_shoot; layout= (1 , 2 ), leftmargin= 15 px, bottommargin= 15 px, kwargs... )
127140end
128141
129142p0_sol
130143
131- plot (S, p0_sol; size= (800 , 450 ))
144+ plot (S, p0_sol; size= (800 , 450 ))
0 commit comments