@@ -84,8 +84,8 @@ savefig("goddard_fig1.png")
84
84
u0 = 0
85
85
u1 = 1
86
86
87
- H0 (x, p) = Lift (F0)(x, p )
88
- H1 (x, p) = Lift (F1)(x, p)
87
+ H0 = Lift (F0); H1 = Lift (F1 )
88
+ # H0(x, p) = Lift(F0)(x, p); H1(x, p) = Lift(F1)(x, p) # debug
89
89
H01 = @Poisson { H0, H1 }
90
90
H001 = @Poisson { H0, H01 }
91
91
H101 = @Poisson { H1, H01 }
@@ -120,16 +120,16 @@ t = direct_sol.times
120
120
x = direct_sol. state
121
121
u = direct_sol. control
122
122
p = direct_sol. costate
123
- H1 (t) = H1 (x (t), p (t))
123
+ φ (t) = H1 (x (t), p (t))
124
124
125
125
u_plot = plot (t, t -> u (t)[1 ], label = " u(t)" )
126
- H1_plot = plot (t, H1 , label = " H₁(x(t), p(t))" )
126
+ H1_plot = plot (t, φ , label = " H₁(x(t), p(t))" )
127
127
g_plot = plot (t, g ∘ x , label = " g(x(t))" )
128
128
display (plot (u_plot, H1_plot, g_plot, layout= (3 ,1 ), size= (700 ,600 )))
129
129
savefig (" goddard_fig2.png" )
130
130
131
131
η = 1e-3
132
- t13 = t[ abs .(H1 .(t)) .≤ η ]
132
+ t13 = t[ abs .(φ .(t)) .≤ η ]
133
133
t23 = t[ 0 .≤ (g ∘ x). (t) .≤ η ]
134
134
p0 = p (t0)
135
135
t1 = min (t13... )
0 commit comments