From 3dcfcab851c5b3330fcc489bc5f41521c945f492 Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Tue, 27 Aug 2024 09:22:14 +0200 Subject: [PATCH] build based on 1777a7d --- dev/.documenter-siteinfo.json | 2 +- ...e-4eab1c87.svg => 2D-example-80f30b14.svg} | 192 +- ...e-110720c7.svg => 2D-example-8181b506.svg} | 192 +- ...e-050c0c36.svg => 2D-example-eacd5bfb.svg} | 192 +- ...e-8e614b4e.svg => 2D-example-fc48b0e7.svg} | 58 +- dev/2D-example.html | 142 +- ...4b7.svg => 2D-preconditioner-2a332038.svg} | 320 +- ...d8b.svg => 2D-preconditioner-867a4fce.svg} | 254 +- ...e71.svg => 2D-preconditioner-d1bc11f9.svg} | 4170 ++++++++--------- ...a41.svg => 2D-preconditioner-e31b3f12.svg} | 328 +- dev/2D-preconditioner.html | 82 +- dev/index.html | 6 +- 12 files changed, 2969 insertions(+), 2969 deletions(-) rename dev/{2D-example-4eab1c87.svg => 2D-example-80f30b14.svg} (84%) rename dev/{2D-example-110720c7.svg => 2D-example-8181b506.svg} (78%) rename dev/{2D-example-050c0c36.svg => 2D-example-eacd5bfb.svg} (78%) rename dev/{2D-example-8e614b4e.svg => 2D-example-fc48b0e7.svg} (86%) rename dev/{2D-preconditioner-a9b004b7.svg => 2D-preconditioner-2a332038.svg} (85%) rename dev/{2D-preconditioner-be377d8b.svg => 2D-preconditioner-867a4fce.svg} (81%) rename dev/{2D-preconditioner-5b03ae71.svg => 2D-preconditioner-d1bc11f9.svg} (64%) rename dev/{2D-preconditioner-2953ea41.svg => 2D-preconditioner-e31b3f12.svg} (85%) diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 9f8b654..0345fcd 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-26T16:49:34","documenter_version":"1.6.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-27T09:22:05","documenter_version":"1.6.0"}} \ No newline at end of file diff --git a/dev/2D-example-4eab1c87.svg b/dev/2D-example-80f30b14.svg similarity index 84% rename from dev/2D-example-4eab1c87.svg rename to dev/2D-example-80f30b14.svg index 26d3744..b57b8c2 100644 --- a/dev/2D-example-4eab1c87.svg +++ b/dev/2D-example-80f30b14.svg @@ -1,118 +1,118 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-110720c7.svg b/dev/2D-example-8181b506.svg similarity index 78% rename from dev/2D-example-110720c7.svg rename to dev/2D-example-8181b506.svg index 1df0682..02647fa 100644 --- a/dev/2D-example-110720c7.svg +++ b/dev/2D-example-8181b506.svg @@ -1,118 +1,118 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-050c0c36.svg b/dev/2D-example-eacd5bfb.svg similarity index 78% rename from dev/2D-example-050c0c36.svg rename to dev/2D-example-eacd5bfb.svg index 864d11e..4c5871d 100644 --- a/dev/2D-example-050c0c36.svg +++ b/dev/2D-example-eacd5bfb.svg @@ -1,118 +1,118 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-8e614b4e.svg b/dev/2D-example-fc48b0e7.svg similarity index 86% rename from dev/2D-example-8e614b4e.svg rename to dev/2D-example-fc48b0e7.svg index 29fbfda..def47e3 100644 --- a/dev/2D-example-8e614b4e.svg +++ b/dev/2D-example-fc48b0e7.svg @@ -1,44 +1,44 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example.html b/dev/2D-example.html index 3096e56..36e776e 100644 --- a/dev/2D-example.html +++ b/dev/2D-example.html @@ -32,45 +32,45 @@ # Plot plot(range(-7, 2, 500), S, xlim = [-7, 2]) plot!([-7,2], [0,0], color = :black) -plot!(xlabel = "p0", ylabel = "S", legend=false)Example block output

Finite difference method

The main goal now is to find the zero of $S$. To this purpose, we use the numerical solver $\texttt{hybrd1}$ given in the package $\texttt{MINPACK.jl}$. If we don't provide the Jacobian $J_S$ of $S$ to the solver, the finite difference method is used to approximate it. 

ξ = [-1.0]                                                  # initial guess
+plot!(xlabel = "p0", ylabel = "S", legend=false)
Example block output

Finite difference method

The main goal now is to find the zero of $S$. To this purpose, we use the numerical solver $\texttt{hybrd1}$ given in the package $\texttt{MINPACK.jl}$. If we don't provide the Jacobian $J_S$ of $S$ to the solver, the finite difference method is used to approximate it. 

ξ = [-1.0]                                                  # initial guess
 S!(s, ξ) = (s[:] .= S(ξ[1]); nothing)                       # intermediate function
 p0_sol = fsolve(S!, ξ, show_trace = true)                   # solve
 println(p0_sol)
Iter     f(x) inf-norm    Step 2-norm      Step time
 ------   --------------   --------------   --------------
-     1     3.000000e+00     0.000000e+00         0.169750
-     2     9.318967e-01     1.069263e+00         0.231582
-     3     3.908179e-08     2.171079e-01         0.001024
-     4     3.126006e-08     3.818465e-16         0.000934
-     5     2.602624e-08     7.541206e-17         0.000891
-     6     1.122729e-07     1.864781e-15         0.000932
-     7     1.251064e-08     1.228963e-15         0.000933
-     8     1.071865e-08     6.960129e-18         0.000945
-     9     3.372047e-08     2.490162e-16         0.001094
-    10     6.149611e-09     1.433785e-16         0.000929
-    11     2.836331e-08     2.624355e-17         0.000876
-    12     2.550366e-08     1.772446e-17         0.000883
-    13     2.794534e-08     1.292187e-17         0.001796
-    14     2.543587e-08     1.364932e-17         0.000886
-    15     6.042093e-09     4.286522e-19         0.001122
-    16     3.036848e-08     7.912620e-17         0.000871
-    17     2.599380e-08     5.504418e-17         0.000881
-    18     2.518026e-08     1.434518e-18         0.000865
-    19     6.005500e-09     5.041160e-20         0.000937
-    20     3.036868e-08     7.817582e-17         0.000882
-    21     2.601931e-08     5.449264e-17         0.000882
-    22     2.521370e-08     1.406704e-18         0.001080
-    23     5.969729e-09     4.888090e-20         0.000940
-    24     3.036827e-08     7.723669e-17         0.000886
-    25     2.604422e-08     5.394382e-17         0.000904
-    26     2.524639e-08     1.379606e-18         0.000862
-    27     5.934758e-09     4.741150e-20         0.001094
-    28     3.036860e-08     7.634290e-17         0.000882
-    29     2.606871e-08     5.342255e-17         0.000877
-    30     2.527841e-08     1.353711e-18         0.001016
-    31     5.900549e-09     4.601691e-20         0.000930
-    32     3.036849e-08     7.546245e-17         0.000860
-    33     2.609265e-08     5.290603e-17         0.000885
-    34     2.530975e-08     1.328510e-18         0.001026
+     1     3.000000e+00     0.000000e+00         0.338011
+     2     9.318967e-01     1.069263e+00         0.351585
+     3     3.908179e-08     2.171079e-01         0.001855
+     4     3.126006e-08     3.818465e-16         0.001820
+     5     2.602624e-08     7.541206e-17         0.001713
+     6     1.122729e-07     1.864781e-15         0.001676
+     7     1.251064e-08     1.228963e-15         0.001689
+     8     1.071865e-08     6.960129e-18         0.001540
+     9     3.372047e-08     2.490162e-16         0.001534
+    10     6.149611e-09     1.433785e-16         0.001656
+    11     2.836331e-08     2.624355e-17         0.001569
+    12     2.550366e-08     1.772446e-17         0.001811
+    13     2.794534e-08     1.292187e-17         0.003071
+    14     2.543587e-08     1.364932e-17         0.001489
+    15     6.042093e-09     4.286522e-19         0.001773
+    16     3.036848e-08     7.912620e-17         0.001765
+    17     2.599380e-08     5.504418e-17         0.001707
+    18     2.518026e-08     1.434518e-18         0.001598
+    19     6.005500e-09     5.041160e-20         0.001717
+    20     3.036868e-08     7.817582e-17         0.001433
+    21     2.601931e-08     5.449264e-17         0.001588
+    22     2.521370e-08     1.406704e-18         0.001452
+    23     5.969729e-09     4.888090e-20         0.001601
+    24     3.036827e-08     7.723669e-17         0.001547
+    25     2.604422e-08     5.394382e-17         0.001487
+    26     2.524639e-08     1.379606e-18         0.001563
+    27     5.934758e-09     4.741150e-20         0.001551
+    28     3.036860e-08     7.634290e-17         0.001513
+    29     2.606871e-08     5.342255e-17         0.001476
+    30     2.527841e-08     1.353711e-18         0.001654
+    31     5.900549e-09     4.601691e-20         0.001646
+    32     3.036849e-08     7.546245e-17         0.001450
+    33     2.609265e-08     5.290603e-17         0.001513
+    34     2.530975e-08     1.328510e-18         0.001515
 Results of Nonlinear Solver Algorithm
  * Algorithm: Modified Powell
  * Starting Point: [-1.0]
@@ -78,10 +78,10 @@
  * Inf-norm of residuals: 0.000000
  * Convergence: true
  * Message: algorithm estimates that the relative error between x and the solution is at most tol
- * Total time: 0.432163 seconds
+ * Total time: 0.742641 seconds
  * Function Calls: 34
  * Jacobian Calls (df/dx): 2
sol = ϕ((t0, tf), x0, p0_sol.x)                             # get the optimal trajectory
-plot(sol)                                                   # plot
Example block output

Automatic differentiation (wrong way)

Now, we want to provide $J_S$ to the solver, thanks to the $\texttt{ForwardDiff.jl}$ package. This Jacobian is computed with the variational equation, and leads to a false result in our case.

+plot(sol) # plotExample block output

Automatic differentiation (wrong way)

Now, we want to provide $J_S$ to the solver, thanks to the $\texttt{ForwardDiff.jl}$ package. This Jacobian is computed with the variational equation, and leads to a false result in our case.

Details. @@ -138,17 +138,17 @@ p0_sol = fsolve(S!, JS!, ξ, show_trace = true) # solve println(p0_sol)
Iter     f(x) inf-norm    Step 2-norm      Step time
 ------   --------------   --------------   --------------
-     1     3.000000e+00     0.000000e+00         0.006896
-     2     5.000000e+00     1.000000e+04         0.064947
-     3     5.000000e+00     3.906250e+03         0.000169
-     4     5.000000e+00     3.906250e+03         0.001459
-     5     5.000000e+00     3.906250e+03         0.000102
-     6     5.000000e+00     5.493164e+02         0.000083
-     7     5.000000e+00     7.724762e+01         0.000107
-     8     9.550781e-01     1.086295e+01         0.001146
-     9     2.633355e-08     2.280435e-01         0.000965
-    10     4.034347e-08     1.733639e-16         0.001006
-    11     7.549755e-09     6.346772e-17         0.000997
+     1     3.000000e+00     0.000000e+00         0.010767
+     2     5.000000e+00     1.000000e+04         0.093134
+     3     5.000000e+00     3.906250e+03         0.000450
+     4     5.000000e+00     3.906250e+03         0.002625
+     5     5.000000e+00     3.906250e+03         0.000278
+     6     5.000000e+00     5.493164e+02         0.000159
+     7     5.000000e+00     7.724762e+01         0.000303
+     8     9.550781e-01     1.086295e+01         0.002040
+     9     2.633355e-08     2.280435e-01         0.001879
+    10     4.034347e-08     1.733639e-16         0.001972
+    11     7.549755e-09     6.346772e-17         0.001704
 Results of Nonlinear Solver Algorithm
  * Algorithm: Modified Powell (User Jac, Expert)
  * Starting Point: [-1.0]
@@ -156,10 +156,10 @@
  * Inf-norm of residuals: 0.000000
  * Convergence: true
  * Message: algorithm estimates that the relative error between x and the solution is at most tol
- * Total time: 0.077891 seconds
+ * Total time: 0.115346 seconds
  * Function Calls: 11
  * Jacobian Calls (df/dx): 2
sol = ϕ((t0, tf), x0, p0_sol.x)                             # get the optimal trajectory
-plt = plot(sol)                                             # plot
Example block output

Automatic differentiation (good way)

The goal is to provide the true Jacobian of $S$ by using the $\texttt{ForwardDiff}$ package, and so we need to indicate to the solver that the dynamic of the system change when $p = 0$.

To understand why we need to give this information to the solver, see the following details.

+plt = plot(sol) # plotExample block output

Automatic differentiation (good way)

The goal is to provide the true Jacobian of $S$ by using the $\texttt{ForwardDiff}$ package, and so we need to indicate to the solver that the dynamic of the system change when $p = 0$.

To understand why we need to give this information to the solver, see the following details.

Details. @@ -245,25 +245,25 @@ JS(ξ) : 2.0000000000000018 Iter f(x) inf-norm Step 2-norm Step time ------ -------------- -------------- -------------- - 1 3.000000e+00 0.000000e+00 1.367424 - 2 6.994405e-15 2.250000e+00 0.008755 - 3 4.551914e-15 1.262177e-29 0.000217 - 4 3.441691e-15 1.774937e-30 0.000147 - 5 2.775558e-15 1.972152e-31 0.000144 - 6 3.441691e-15 1.972152e-31 0.000132 - 7 2.775558e-15 1.972152e-31 0.000138 - 8 2.997602e-15 1.774937e-30 0.000318 - 9 1.443290e-15 1.972152e-31 0.000167 - 10 2.997602e-15 7.888609e-31 0.000137 - 11 2.997602e-15 1.972152e-31 0.000135 - 12 2.997602e-15 1.972152e-31 0.000268 - 13 2.997602e-15 1.972152e-31 0.000140 - 14 1.443290e-15 1.972152e-31 0.000129 - 15 2.343750e-02 1.373291e-04 0.000146 - 16 2.997602e-15 1.373291e-04 0.000131 - 17 2.997602e-15 1.972152e-31 0.000137 - 18 1.443290e-15 1.972152e-31 0.000130 - 19 1.464844e-03 5.364418e-07 0.000135 + 1 3.000000e+00 0.000000e+00 2.068230 + 2 6.994405e-15 2.250000e+00 0.013056 + 3 4.551914e-15 1.262177e-29 0.000390 + 4 3.441691e-15 1.774937e-30 0.000327 + 5 2.775558e-15 1.972152e-31 0.000320 + 6 3.441691e-15 1.972152e-31 0.000298 + 7 2.775558e-15 1.972152e-31 0.000308 + 8 2.997602e-15 1.774937e-30 0.000634 + 9 1.443290e-15 1.972152e-31 0.088895 + 10 2.997602e-15 7.888609e-31 0.000699 + 11 2.997602e-15 1.972152e-31 0.000341 + 12 2.997602e-15 1.972152e-31 0.000793 + 13 2.997602e-15 1.972152e-31 0.000367 + 14 1.443290e-15 1.972152e-31 0.000339 + 15 2.343750e-02 1.373291e-04 0.000334 + 16 2.997602e-15 1.373291e-04 0.000271 + 17 2.997602e-15 1.972152e-31 0.000318 + 18 1.443290e-15 1.972152e-31 0.000262 + 19 1.464844e-03 5.364418e-07 0.000313 Results of Nonlinear Solver Algorithm * Algorithm: Modified Powell (User Jac, Expert) * Starting Point: [-1.0] @@ -271,7 +271,7 @@ * Inf-norm of residuals: 0.000000 * Convergence: true * Message: iteration is not making good progress, measured by improvement from last 10 iterations - * Total time: 1.378944 seconds + * Total time: 2.176545 seconds * Function Calls: 19 * Jacobian Calls (df/dx): 3
# get optimal trajectory
 sol_ = φ((t0, tf), x0, p0_sol.x[1], saveat=range(t0, tf, 500))
@@ -288,4 +288,4 @@
 plt_u = plot(t, u, label = "u")
 
 plt_xp = plot(plt_x, plt_p, layout=(1, 2))
-plot(plt_xp, plt_u, layout = (2, 1))
Example block output
+plot(plt_xp, plt_u, layout = (2, 1))Example block output
diff --git a/dev/2D-preconditioner-a9b004b7.svg b/dev/2D-preconditioner-2a332038.svg similarity index 85% rename from dev/2D-preconditioner-a9b004b7.svg rename to dev/2D-preconditioner-2a332038.svg index 64284c0..789881e 100644 --- a/dev/2D-preconditioner-a9b004b7.svg +++ b/dev/2D-preconditioner-2a332038.svg @@ -1,190 +1,190 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-preconditioner-be377d8b.svg b/dev/2D-preconditioner-867a4fce.svg similarity index 81% rename from dev/2D-preconditioner-be377d8b.svg rename to dev/2D-preconditioner-867a4fce.svg index 919e4f5..8e8e4db 100644 --- a/dev/2D-preconditioner-be377d8b.svg +++ b/dev/2D-preconditioner-867a4fce.svg @@ -1,146 +1,146 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-preconditioner-5b03ae71.svg b/dev/2D-preconditioner-d1bc11f9.svg similarity index 64% rename from dev/2D-preconditioner-5b03ae71.svg rename to dev/2D-preconditioner-d1bc11f9.svg index c7f43b2..21e7261 100644 --- a/dev/2D-preconditioner-5b03ae71.svg +++ b/dev/2D-preconditioner-d1bc11f9.svg @@ -1,2107 +1,2107 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-preconditioner-2953ea41.svg b/dev/2D-preconditioner-e31b3f12.svg similarity index 85% rename from dev/2D-preconditioner-2953ea41.svg rename to dev/2D-preconditioner-e31b3f12.svg index b0e6eca..4daf6c5 100644 --- a/dev/2D-preconditioner-2953ea41.svg +++ b/dev/2D-preconditioner-e31b3f12.svg @@ -1,194 +1,194 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-preconditioner.html b/dev/2D-preconditioner.html index 7706808..1f8a695 100644 --- a/dev/2D-preconditioner.html +++ b/dev/2D-preconditioner.html @@ -93,23 +93,23 @@ p0_sol = fsolve(S₂!, JS₂!, ξ, show_trace = true) # solve println(p0_sol) # print solution
Iter     f(x) inf-norm    Step 2-norm      Step time
 ------   --------------   --------------   --------------
-     1     3.845299e+00     0.000000e+00         0.003686
-     2     5.000000e+00     1.559496e+00         5.050275
-     3     5.000000e+00     4.983094e-01         0.000091
-     4     5.000000e+00     4.983094e-01         0.000387
-     5     5.000000e+00     4.983094e-01         0.000034
-     6     2.825589e+00     9.417471e-02         0.000296
-     7     3.275582e+00     5.569993e-02         0.000281
-     8     1.835694e+00     1.605481e-02         0.000216
-     9     5.000000e+00     4.108358e-02         0.000039
-    10     8.767661e-01     2.198078e-02         0.000180
-    11     1.052824e+00     2.476848e-03         0.000227
-    12     2.252319e-01     7.373618e-04         0.000179
-    13     6.733089e-02     6.111140e-05         0.000255
-    14     4.061534e-03     3.236777e-06         0.000175
-    15     6.974839e-05     1.047584e-08         0.000244
-    16     7.333033e-08     3.198326e-12         0.000171
-    17     1.322831e-12     3.527840e-18         0.000246
+     1     3.845299e+00     0.000000e+00         0.005766
+     2     5.000000e+00     1.559496e+00         8.658186
+     3     5.000000e+00     4.983094e-01         0.000386
+     4     5.000000e+00     4.983094e-01         0.000678
+     5     5.000000e+00     4.983094e-01         0.000152
+     6     2.825589e+00     9.417471e-02         0.000485
+     7     3.275582e+00     5.569993e-02         0.000410
+     8     1.835694e+00     1.605481e-02         0.000444
+     9     5.000000e+00     4.108358e-02         0.000075
+    10     8.767661e-01     2.198078e-02         0.000323
+    11     1.052824e+00     2.476848e-03         0.000341
+    12     2.252319e-01     7.373618e-04         0.000328
+    13     6.733089e-02     6.111140e-05         0.000378
+    14     4.061534e-03     3.236777e-06         0.000309
+    15     6.974839e-05     1.047584e-08         0.000349
+    16     7.333033e-08     3.198326e-12         0.000302
+    17     1.322831e-12     3.527840e-18         0.000356
 Results of Nonlinear Solver Algorithm
  * Algorithm: Modified Powell (User Jac, Expert)
  * Starting Point: [-0.5]
@@ -117,7 +117,7 @@
  * Inf-norm of residuals: 0.000000
  * Convergence: true
  * Message: algorithm estimates that the relative error between x and the solution is at most tol
- * Total time: 5.057001 seconds
+ * Total time: 8.669299 seconds
  * Function Calls: 17
  * Jacobian Calls (df/dx): 2
sol = ϕ((t0, tf), x0, [η(p0_sol.x[1]), p0_sol.x[1]],
     saveat=range(t0, tf, 500))                                  # get optimal trajectory
@@ -137,7 +137,7 @@
 plt_u  = plot(t, u,  label = "u" )
 
 plt_xp = plot(plt_x⁰, plt_p⁰, plt_x, plt_p, layout=(2, 2))
-plot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))
Example block output

Construction of the geometric preconditioner

The goal is now to use the geometric preconditioning method proposed in [mettre article]. For this purpose, the first step is to create some points in the boundary of the accessible augmented set, and to fit an ellipse on these points.

The second step is to create the linear diffeomorphism $\phi \colon \mathbb R^2 \to \mathbb R^2, (x^0,x_1) \to Ax + B$ that transforms the fitted ellipse into the unit circle, and that satisfy the condition

\[ \frac{\partial \phi}{\partial x^0} = k e_1, \]

with $k>0$ and $e_1 = (1,0)$. In this context, we denote

\[ A = \left( \begin{array}{cc} k & A_{x^0} \\ 0 & A_{x_1} \end{array} \right) +plot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))Example block output

Construction of the geometric preconditioner

The goal is now to use the geometric preconditioning method proposed in [mettre article]. For this purpose, the first step is to create some points in the boundary of the accessible augmented set, and to fit an ellipse on these points.

The second step is to create the linear diffeomorphism $\phi \colon \mathbb R^2 \to \mathbb R^2, (x^0,x_1) \to Ax + B$ that transforms the fitted ellipse into the unit circle, and that satisfy the condition

\[ \frac{\partial \phi}{\partial x^0} = k e_1, \]

with $k>0$ and $e_1 = (1,0)$. In this context, we denote

\[ A = \left( \begin{array}{cc} k & A_{x^0} \\ 0 & A_{x_1} \end{array} \right) \qquad \text{and} \qquad B = \left( \begin{array}{c} B_{x^0} \\ B_{x_1} \end{array} \right).\]

This diffeomorphism is given from the semi-axis $a,b >0$, the angle $\theta \in [0, \frac{\pi}{2}[$ between the semi-axis $b$ and the $x$-axis, and the center $c \in \mathbb R^2$ by

\[ \phi(x) = r(-\beta_0) s(a^{-1}, b^{-1}) r(\theta) (x - c),\]

where $r$ and $s$ correspond respectively to the rotation and the scale matrix, defined by

\[ r(\theta) = \left( \begin{array}{cc} \phantom - \cos(\theta) & \sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{array} \right) \qquad \text{and} \qquad @@ -202,7 +202,7 @@ plot!(yₑ[2,:], yₑ[1,:], label = "") plot!( xlabel = "y", ylabel = "y⁰") -plot(plt_x, plt_y, layout = (1,2), size=(800, 400))Example block output

The general shooting function $T \colon \mathbb R^2 \to \mathbb R$ in the new coordinate is defined by

\[ T(q) = A_{x_1} \varphi \big(t_0, x_0, p_0(A^\top q), t_f) + B_{x_1} - y_T,\]

where the function $p_0 \colon \mathbb R^2 \to \mathbb R^2$ correspond to the mapping between the final and the initial augmented costate, and is given by

\[ p_0(p^0, p_1) = (p^0, p_1 + p^0 t_f),\]

and $y_T = A_{x_1} x_T + B_x$ to the target in the new system of coordinates. By using the definition of $S$ and $y_T$, we obtain

\[\begin{align*} +plot(plt_x, plt_y, layout = (1,2), size=(800, 400))Example block output

The general shooting function $T \colon \mathbb R^2 \to \mathbb R$ in the new coordinate is defined by

\[ T(q) = A_{x_1} \varphi \big(t_0, x_0, p_0(A^\top q), t_f) + B_{x_1} - y_T,\]

where the function $p_0 \colon \mathbb R^2 \to \mathbb R^2$ correspond to the mapping between the final and the initial augmented costate, and is given by

\[ p_0(p^0, p_1) = (p^0, p_1 + p^0 t_f),\]

and $y_T = A_{x_1} x_T + B_x$ to the target in the new system of coordinates. By using the definition of $S$ and $y_T$, we obtain

\[\begin{align*} T(q) &= A_{x_1} \varphi \big(t_0, x_0, p_0(A^\top q), t_f) + B_{x_1} - (A_{x_1} x_T + B_{x_1}) \\ &= A_{x_1} (S \circ p_0)(A^\top q). \end{align*}\]

which highlight that the proposed geometric preconditioning method is a left and right side preconditioner of the shooting function. Finally, we define the two shooting functions $T_1$ and $T_2$ by using the two method of normalization used before for the function $S$

\[ T_1(p) = T(-1, p) \qquad \text{and} \qquad T_2(p) = T \big(\eta(p), p \big).\]

p₀(p) = [p[1], p[2] + p[1]*tf]                                  # function p₀(⋅)
@@ -248,24 +248,24 @@
 q_sol = fsolve(T₂!, JT₂!, ξ, show_trace = true)                 # solve
 println(q_sol)                                                  # print solution
Iter     f(x) inf-norm    Step 2-norm      Step time
 ------   --------------   --------------   --------------
-     1     3.223430e-01     0.000000e+00         0.000475
-     2     7.034101e-02     1.406250e-01         5.073508
-     3     1.135104e-02     1.095650e-02         0.000320
-     4     1.033639e-04     4.056828e-04         0.000213
-     5     2.155036e-08     3.426091e-08         0.000201
-     6     8.398089e-16     1.489880e-15         0.000182
-     7     3.119290e-16     2.262576e-30         0.000185
-     8     2.399454e-17     1.659615e-31         0.000176
-     9     2.399454e-17     1.152511e-33         0.000183
-    10     3.119290e-16     1.947743e-31         0.000210
-    11     2.399454e-17     2.028246e-31         0.000505
-    12     2.399454e-17     2.097361e-32         0.000260
-    13     2.399454e-17     2.593149e-33         0.000217
-    14     2.399454e-17     6.482873e-34         0.000176
-    15     2.399454e-17     1.620718e-34         0.000179
-    16     2.399454e-17     4.051796e-35         0.000176
-    17     2.399454e-17     1.012949e-35         0.000181
-    18     2.399454e-17     2.532372e-36         0.000170
+     1     3.223430e-01     0.000000e+00         0.001437
+     2     7.034101e-02     1.406250e-01         8.213185
+     3     1.135104e-02     1.095650e-02         0.000532
+     4     1.033639e-04     4.056828e-04         0.000458
+     5     2.155036e-08     3.426091e-08         0.000381
+     6     8.398089e-16     1.489880e-15         0.000353
+     7     3.119290e-16     2.262576e-30         0.000388
+     8     2.399454e-17     1.659615e-31         0.000342
+     9     2.399454e-17     1.152511e-33         0.000359
+    10     3.119290e-16     1.947743e-31         0.000350
+    11     2.399454e-17     2.028246e-31         0.000964
+    12     2.399454e-17     2.097361e-32         0.000348
+    13     2.399454e-17     2.593149e-33         0.000418
+    14     2.399454e-17     6.482873e-34         0.000371
+    15     2.399454e-17     1.620718e-34         0.000375
+    16     2.399454e-17     4.051796e-35         0.000411
+    17     2.399454e-17     1.012949e-35         0.000356
+    18     2.399454e-17     2.532372e-36         0.000344
 Results of Nonlinear Solver Algorithm
  * Algorithm: Modified Powell (User Jac, Expert)
  * Starting Point: [0.5]
@@ -273,7 +273,7 @@
  * Inf-norm of residuals: 0.000000
  * Convergence: true
  * Message: iteration is not making good progress, measured by improvement from last 10 iterations
- * Total time: 5.077542 seconds
+ * Total time: 8.221452 seconds
  * Function Calls: 18
  * Jacobian Calls (df/dx): 2
p0_sol = p₀(transpose(A)*[η(q_sol.x[1]), q_sol.x[1]])           # get the optimal initial costate in old coordinates
 sol = ϕ((t0, tf), x0, p0_sol, saveat=range(t0, tf, 500))        # get optimal trajectory
@@ -293,7 +293,7 @@
 plt_u  = plot(t, u,  label = "u" )
 
 plt_xp = plot(plt_x⁰, plt_p⁰, plt_x, plt_p, layout=(2, 2))
-plot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))
Example block output

Comparison

It is shown in [mettre article] that if the boundary of augmented accessible set is the fitted ellipse then the shooting function $T_2$ is the identity function. Due to the error of the approximation, the function $T_2$ is not the identity, but we hope that this function is close to this ideal function, and so that the convergence of $T_2$ is faster than the one of $S_2$.

The following code compares the convergence of these two shooting function. We also study for these functions the domain of initial guesses for the solver that hit the bounds $[-1, 1]$ during the solving process. 

## initial guesses
+plot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))
Example block output

Comparison

It is shown in [mettre article] that if the boundary of augmented accessible set is the fitted ellipse then the shooting function $T_2$ is the identity function. Due to the error of the approximation, the function $T_2$ is not the identity, but we hope that this function is close to this ideal function, and so that the convergence of $T_2$ is faster than the one of $S_2$.

The following code compares the convergence of these two shooting function. We also study for these functions the domain of initial guesses for the solver that hit the bounds $[-1, 1]$ during the solving process. 

## initial guesses
 N = 1000; ξ = range(-1.,1.,N)
 
 # initialization of the matrix
@@ -368,4 +368,4 @@
 scatter!(1,1, markercolor = :red, markerstrokecolor = :red, marker = 2, label = "not converged")
 
 plt_2 = plot(plt_21, plt_22, layout = (1,2))
-plot(plt_1, plt_2, layout = grid(2,1, heights = [0.5, 0.5]), size=(700, 600))
Example block output
+plot(plt_1, plt_2, layout = grid(2,1, heights = [0.5, 0.5]), size=(700, 600))Example block output
diff --git a/dev/index.html b/dev/index.html index f9f726a..36293ea 100644 --- a/dev/index.html +++ b/dev/index.html @@ -5,11 +5,11 @@ \phantom{\mathrm{s.c.}~} u(t) \in [-1,1], & t\in [t_0, t_f], \\[0.5em] \phantom{\mathrm{s.c.}~} x(t_0) = x_0, \quad x(t_f) = x_f, \end{array} \right.\]

with $x_0$, $t_0$, $x_f$ and $t_f$ fixed. This problem is simple, and can be analytically solve without the use of numerical method. However, the goal is to solve this problem by indirect shooting.

Dependencies

All the numerical simulations to generate this documentation are performed with the following packages.

using Pkg
-Pkg.status()
Status `~/work/preconditioning/preconditioning/docs/Project.toml`
+Pkg.status()
Status `~/actions-runner/_work/preconditioning/preconditioning/docs/Project.toml`
   [0c46a032] DifferentialEquations v7.13.0
   [e30172f5] Documenter v1.6.0
   [f6369f11] ForwardDiff v0.10.36
-  [4fb4c77a] GeometricPreconditioner v0.1.0 `~/work/preconditioning/preconditioning`
+  [4fb4c77a] GeometricPreconditioner v0.1.0 `~/actions-runner/_work/preconditioning/preconditioning`
   [4854310b] MINPACK v1.3.0
   [5f98b655] OptimalControl v0.11.4
-  [91a5bcdd] Plots v1.40.5
+ [91a5bcdd] Plots v1.40.7