From d08724a4651c0f9a8ff1af4326dd8e9997ece574 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 29 Aug 2023 17:37:01 +0900 Subject: [PATCH 01/11] Minor fix in docs. --- docs/src/modules/integrators.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/src/modules/integrators.md b/docs/src/modules/integrators.md index ad04e058e..a71e1cf29 100644 --- a/docs/src/modules/integrators.md +++ b/docs/src/modules/integrators.md @@ -78,8 +78,10 @@ Pages = ["integrators/rk/abstract.jl", ```@autodocs Modules = [GeometricIntegrators.Integrators] -Pages = ["integrators/vi/integrators_vprk.jl"] -Pages = ["integrators/vi/vprk_methods.jl"] +Pages = [ + "integrators/vi/integrators_vprk.jl", + "integrators/vi/vprk_methods.jl", + ] ``` ## Degenerate Variational Integrators From b23b13e24531788321638abe35dca7fe97705147 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 29 Aug 2023 17:41:20 +0900 Subject: [PATCH 02/11] Rename VPRK files. --- docs/src/modules/integrators.md | 2 +- src/Integrators.jl | 4 ++-- .../vi/{integrators_vprk_cache.jl => vprk_cache.jl} | 0 .../vi/{integrators_vprk.jl => vprk_integrator.jl} | 0 4 files changed, 3 insertions(+), 3 deletions(-) rename src/integrators/vi/{integrators_vprk_cache.jl => vprk_cache.jl} (100%) rename src/integrators/vi/{integrators_vprk.jl => vprk_integrator.jl} (100%) diff --git a/docs/src/modules/integrators.md b/docs/src/modules/integrators.md index a71e1cf29..f60f5aa37 100644 --- a/docs/src/modules/integrators.md +++ b/docs/src/modules/integrators.md @@ -79,7 +79,7 @@ Pages = ["integrators/rk/abstract.jl", ```@autodocs Modules = [GeometricIntegrators.Integrators] Pages = [ - "integrators/vi/integrators_vprk.jl", + "integrators/vi/vprk_integrator.jl", "integrators/vi/vprk_methods.jl", ] ``` diff --git a/src/Integrators.jl b/src/Integrators.jl index 7b65dabf0..9874174f8 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -192,8 +192,8 @@ module Integrators export IntegratorVPRK include("integrators/vi/vprk_methods.jl") - include("integrators/vi/integrators_vprk_cache.jl") - include("integrators/vi/integrators_vprk.jl") + include("integrators/vi/vprk_cache.jl") + include("integrators/vi/vprk_integrator.jl") # export IntegratorCGVI, IntegratorDGVI, IntegratorDGVIEXP, diff --git a/src/integrators/vi/integrators_vprk_cache.jl b/src/integrators/vi/vprk_cache.jl similarity index 100% rename from src/integrators/vi/integrators_vprk_cache.jl rename to src/integrators/vi/vprk_cache.jl diff --git a/src/integrators/vi/integrators_vprk.jl b/src/integrators/vi/vprk_integrator.jl similarity index 100% rename from src/integrators/vi/integrators_vprk.jl rename to src/integrators/vi/vprk_integrator.jl From 504ed8116b344748e6f8f5f5e3f6da320ade394f Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 29 Aug 2023 18:42:55 +0900 Subject: [PATCH 03/11] Add skeletons for variational and Hamilton-Pontryagin integrators. --- docs/make.jl | 1 + docs/src/index.md | 58 ++++++++++++------- docs/src/integrators/hpi.md | 1 + docs/src/modules/integrators.md | 36 +++++++----- src/Integrators.jl | 10 ++++ src/integrators/hpi/hpi_methods.jl | 2 + src/integrators/hpi/hpi_midpoint.jl | 50 ++++++++++++++++ src/integrators/hpi/hpi_trapezoidal.jl | 51 ++++++++++++++++ src/integrators/method_list.jl | 8 +++ .../vi/position_momentum_midpoint.jl | 28 +++++++++ .../vi/position_momentum_trapezoidal.jl | 28 +++++++++ src/integrators/vi/vi_methods.jl | 5 ++ src/integrators/vi/vprk_methods.jl | 2 +- .../hamilton_pontryagin_integrators_tests.jl | 0 test/runtests.jl | 2 + 15 files changed, 246 insertions(+), 36 deletions(-) create mode 100644 docs/src/integrators/hpi.md create mode 100644 src/integrators/hpi/hpi_methods.jl create mode 100644 src/integrators/hpi/hpi_midpoint.jl create mode 100644 src/integrators/hpi/hpi_trapezoidal.jl create mode 100644 src/integrators/vi/position_momentum_midpoint.jl create mode 100644 src/integrators/vi/position_momentum_trapezoidal.jl create mode 100644 src/integrators/vi/vi_methods.jl create mode 100644 test/integrators/hamilton_pontryagin_integrators_tests.jl diff --git a/docs/make.jl b/docs/make.jl index 56c71b342..3a51e7fe3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -55,6 +55,7 @@ makedocs(bib, "DVI" => "integrators/dvi.md", # "CGVI" => "integrators/cgvi.md", # "DGVI" => "integrators/dgvi.md", + # "HPI" => "integrators/hpi.md", # "HPG" => "integrators/hpg.md", ], "Modules" => [ diff --git a/docs/src/index.md b/docs/src/index.md index e1e92230a..30416b4e3 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,4 +1,3 @@ - # GeometricIntegrators.jl *Julia library of geometric integrators for differential equations.* @@ -47,38 +46,57 @@ While DifferentialEquations.jl provides a feature-rich ecosystem for the solutio ## Manual ```@contents -Pages = ["tutorial/tutorial.md", - "problems.md", - "methods.md", +Pages = [ + "tutorial/tutorial.md", + "problems.md", + "methods.md", ] ``` -## Modules +## Integrators ```@contents -Pages = ["modules/methods.md", - "modules/integrators.md", - #"modules/discontinuities.md", - "modules/problems.md", - "modules/equations.md", - #"modules/simulations.md", - "modules/solutions.md", - "modules/rungekutta.md", - "modules/rungekutta_partitioned.md", - "modules/spark.md", +Pages = [ + "integrators/rk.md", + "integrators/splitting.md", + "integrators/variational.md", + "integrators/vprk.md", + "integrators/spark.md", + "integrators/dvi.md", + "integrators/cgvi.md", + "integrators/dgvi.md", + "integrators/hpi.md", + "integrators/hpg.md", ] ``` +## Modules + +```@contents +Pages = [ + "modules/methods.md", + "modules/integrators.md", + # "modules/discontinuities.md", + "modules/problems.md", + "modules/equations.md", + # "modules/simulations.md", + "modules/solutions.md", + "modules/rungekutta.md", + "modules/rungekutta_partitioned.md", + "modules/spark.md", +] +``` ## Developer Documentation ```@contents -Pages = ["developer/integrators.md", - "developer/projections.md", - "developer/code_integration.md", - "developer/custom_integrators.md", - "developer/adaptive_time_stepping.md", +Pages = [ + "developer/integrators.md", + "developer/projections.md", + "developer/code_integration.md", + "developer/custom_integrators.md", + "developer/adaptive_time_stepping.md", ] ``` diff --git a/docs/src/integrators/hpi.md b/docs/src/integrators/hpi.md new file mode 100644 index 000000000..7a2b6074b --- /dev/null +++ b/docs/src/integrators/hpi.md @@ -0,0 +1 @@ +# [Hamilton-Pontryagin Integrators](@id hpi) diff --git a/docs/src/modules/integrators.md b/docs/src/modules/integrators.md index f60f5aa37..5df5897d7 100644 --- a/docs/src/modules/integrators.md +++ b/docs/src/modules/integrators.md @@ -2,19 +2,16 @@ # Integrators -## Common +## Geometric Integrator ```@autodocs Modules = [GeometricIntegrators.Integrators] -Pages = ["integrators/abstract_coefficients.jl", - "integrators/abstract_integrator.jl", - "integrators/abstract_tableau.jl", - "integrators/integrator_cache.jl", - "integrators/integrators_common.jl", - "integrators/integrators.jl"] +Pages = [ + "integrators/integrator_cache.jl", + "integrators/integrator.jl", + ] ``` - ## Initial Guesses ```@autodocs @@ -26,7 +23,6 @@ Pages = [ ] ``` - ## Extrapolation Methods The extrapolation routines are exclusively used for computing @@ -42,7 +38,6 @@ Pages = ["integrators/extrapolation/extrapolation.jl", "integrators/extrapolation/midpoint.jl"] ``` - # Euler Integrators ```@autodocs @@ -53,8 +48,7 @@ Pages = [ ] ``` - -## Runge-Kutta Methods +## Runge-Kutta Integrators ```@autodocs Modules = [GeometricIntegrators.Integrators] @@ -73,12 +67,14 @@ Pages = ["integrators/rk/abstract.jl", ] ``` - ## Variational Integrators ```@autodocs Modules = [GeometricIntegrators.Integrators] Pages = [ + "integrators/vi/vi_methods.jl", + "integrators/vi/position_momentum_midpoint.jl", + "integrators/vi/position_momentum_trapezoidal.jl", "integrators/vi/vprk_integrator.jl", "integrators/vi/vprk_methods.jl", ] @@ -89,12 +85,12 @@ Pages = [ ```@autodocs Modules = [GeometricIntegrators.Integrators] Pages = [ + "integrators/dvi/methods.jl", "integrators/dvi/integrators_dvi_a.jl", "integrators/dvi/integrators_dvi_b.jl", "integrators/dvi/integrators_cmdvi.jl", "integrators/dvi/integrators_ctdvi.jl", "integrators/dvi/integrators_dvrk.jl", - "integrators/dvi/methods.jl", ] ``` @@ -106,8 +102,18 @@ Pages = ["integrators/cgvi/integrators_cgvi.jl", "integrators/dgvi/integrators_dgvi.jl"] ``` +## Hamilton-Pontryagin Integrators + +```@autodocs +Modules = [GeometricIntegrators.Integrators] +Pages = [ + "integrators/hpi/hpi_methods.jl", + "integrators/hpi/hpi_midpoint.jl", + "integrators/hpi/hpi_trapezoidal.jl", + ] +``` -## Splitting Methods +## Splitting and Composition Methods ```@autodocs Modules = [GeometricIntegrators.Integrators] diff --git a/src/Integrators.jl b/src/Integrators.jl index 9874174f8..2ffac7bd2 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -191,6 +191,9 @@ module Integrators export IntegratorVPRK + include("integrators/vi/vi_methods.jl") + include("integrators/vi/position_momentum_midpoint.jl") + include("integrators/vi/position_momentum_trapezoidal.jl") include("integrators/vi/vprk_methods.jl") include("integrators/vi/vprk_cache.jl") include("integrators/vi/vprk_integrator.jl") @@ -219,6 +222,13 @@ module Integrators include("integrators/dvi/integrators_dvrk.jl") + export HPImidpointIntegrator, HPItrapezoidalIntegrator + + include("integrators/hpi/hpi_methods.jl") + include("integrators/hpi/hpi_midpoint.jl") + include("integrators/hpi/hpi_trapezoidal.jl") + + export NoProjection, projection include("projections/methods.jl") diff --git a/src/integrators/hpi/hpi_methods.jl b/src/integrators/hpi/hpi_methods.jl new file mode 100644 index 000000000..bd717e6fd --- /dev/null +++ b/src/integrators/hpi/hpi_methods.jl @@ -0,0 +1,2 @@ + +abstract type HPIMethod <: LODEMethod end diff --git a/src/integrators/hpi/hpi_midpoint.jl b/src/integrators/hpi/hpi_midpoint.jl new file mode 100644 index 000000000..a6da90c7c --- /dev/null +++ b/src/integrators/hpi/hpi_midpoint.jl @@ -0,0 +1,50 @@ +@doc raw""" +Hamilton-Pontryagin Integrator using midpoint quadrature. + +We consider a discrete Lagrangian of the form +```math +L_d (q_{n}, q_{n+1}) = h \, L \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) , +``` +where $q_{n}$ approximates the solution $q(t_n)$ and $v_{n+1/2}$ is the velocity, which is assumed to be constant in the interval $[t_{n}, t_{n+1}]$. +The discrete Hamilton-Pontryagin action reads +```math +A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_{d}^{\alpha,\beta} (q_{n}, q_{n+1}, v_{n+1/2}) + \left< p_{n+1/2} , v_{n+1/2} - \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \right> \bigg] , +``` +where $\phi_h$ is a map that computes the velocity $v_{n+1/2}$ as a function of $q_{n}$, $q_{n+1}$ and a set of parameters $a_{n,n+1}$. +A trivial example of such a map that does not depend on any parameters $a_{n,n+1}$ is +```math +\phi_h (q_{n}, q_{n+1}; a_{n,n+1}) = \frac{q_{n+1} - q_{n}}{h} . +``` +In order to solve the discrete Euler-Lagrange equations, the user needs to specify the map $\phi_h$, its gradients with respect to $q_{n}$ and $q_{n+1}$, denoted by $D_1 \phi_h$ and $D_2 \phi_h$, respectively, the gradient with respect to the parameters, denoted by $D_a \phi_h$, and an initial set of parameters $a_{0}$. + +The equations of motion, that are solved by this integrator, is computed as: +```math +\begin{aligned} +0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n-1} + q_{n}}{2}, v_{n-1/2} \bigg) + + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) \\ + &- h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + - h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} , \\ +0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} , \\ +p_{n+1/2} &= \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) , \\ +v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) . +\end{aligned} +``` +The current implementation requires that $D_2 \phi_h $ and $D_a \phi_h$ do not depend on the second and third argument, but only $D_1 \phi_h$ is allowed to depend on all three arguments. +Otherwise it would be necessary to prescribe initial conditions $(q_{-1}, a_{-1,0}, q_{0}, p_{0})$ instead of just $(q_{0}, p_{0})$. + +""" +struct HPImidpoint{ϕT, D₁ϕT, D₂ϕT, DₐϕT, PT} <: HPIMethod + ϕ::ϕT + D₁ϕ::D₁ϕT + D₂ϕ::D₂ϕT + Dₐϕ::DₐϕT + params::PT +end + +isexplicit(method::HPImidpoint) = false +isimplicit(method::HPImidpoint) = true +issymmetric(method::HPImidpoint) = missing +issymplectic(method::HPImidpoint) = true + + +const HPImidpointIntegrator{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:HPImidpoint} diff --git a/src/integrators/hpi/hpi_trapezoidal.jl b/src/integrators/hpi/hpi_trapezoidal.jl new file mode 100644 index 000000000..63dbe226d --- /dev/null +++ b/src/integrators/hpi/hpi_trapezoidal.jl @@ -0,0 +1,51 @@ +@doc raw""" +Hamilton-Pontryagin Integrator using trapezoidal quadrature. + +We consider a discrete Lagrangian of the form +```math +L_d (q_{n}, q_{n+1}) = \frac{h}{2} \big[ L (q_{n}, v_{n+1/2}) + L (q_{n+1}, v_{n+1/2}) \big] , +``` +where $q_{n}$ approximates the solution $q(t_n)$ and $v_{n+1/2}$ is the velocity, which is assumed to be constant in the interval $[t_{n}, t_{n+1}]$. +The discrete Hamilton-Pontryagin action reads +```math +A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_d (q_{n}, q_{n+1}) + h \left< p_{n+1} , v_{n+1/2} - \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \right> \bigg] , +``` +where $\phi_h$ is a map that computes the velocity $v_{n+1/2}$ as a function of $q_{n}$, $q_{n+1}$ and a set of parameters $a_{n,n+1}$. +A trivial example of such a map that does not depend on any parameters $a_{n,n+1}$ is +```math +\phi_h (q_{n}, q_{n+1}; a_{n,n+1}) = \frac{q_{n+1} - q_{n}}{h} . +``` +In order to solve the discrete Euler-Lagrange equations, the user needs to specify the map $\phi_h$, its gradients with respect to $q_{n}$ and $q_{n+1}$, denoted by $D_1 \phi_h$ and $D_2 \phi_h$, respectively, the gradient with respect to the parameters, denoted by $D_a \phi_h$, and an initial set of parameters $a_{0}$. + +The equations of motion, that are solved by this integrator, is computed as: +```math +\begin{aligned} +0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2}) + + \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n+1}, v_{n+1/2}) \\ + &- h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} + - h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n} , \\ +0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} , \\ +p_{n+1} &= \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n}, v_{n+1/2}) + + \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n+1}, v_{n+1/2}) , \\ +v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) . +\end{aligned} +``` +The current implementation requires that $D_2 \phi_h $ and $D_a \phi_h$ do not depend on the second and third argument, but only $D_1 \phi_h$ is allowed to depend on all three arguments. +Otherwise it would be necessary to prescribe initial conditions $(q_{-1}, a_{-1,0}, q_{0}, p_{0})$ instead of just $(q_{0}, p_{0})$. + +""" +struct HPItrapezoidal{ϕT, D₁ϕT, D₂ϕT, DₐϕT, PT} <: HPIMethod + ϕ::ϕT + D₁ϕ::D₁ϕT + D₂ϕ::D₂ϕT + Dₐϕ::DₐϕT + params::PT +end + +isexplicit(method::HPItrapezoidal) = false +isimplicit(method::HPItrapezoidal) = true +issymmetric(method::HPItrapezoidal) = missing +issymplectic(method::HPItrapezoidal) = true + + +const HPItrapezoidalIntegrator{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:HPItrapezoidal} diff --git a/src/integrators/method_list.jl b/src/integrators/method_list.jl index 66f99b2de..70cdd192a 100644 --- a/src/integrators/method_list.jl +++ b/src/integrators/method_list.jl @@ -10,6 +10,8 @@ meta_methods = ( DVRK, ProjectedVPRK, ProjectedMethod, + HPImidpoint, + HPItrapezoidal, ) explicit_rungekutta_methods = ( @@ -124,6 +126,11 @@ variational_partitioned_runge_kutta_families = ( VPRKLobattoIIIGIIIḠ, ) +variational_integrators = ( + PMVImidpoint, + PMVItrapezoidal, +) + degenerate_variational_integrators = ( DVIA, DVIB, @@ -151,6 +158,7 @@ method_groups = ( partitioned_runge_kutta_methods, partitioned_runge_kutta_families, variational_partitioned_runge_kutta_families, + variational_integrators, degenerate_variational_integrators, splitting_methods, ) diff --git a/src/integrators/vi/position_momentum_midpoint.jl b/src/integrators/vi/position_momentum_midpoint.jl new file mode 100644 index 000000000..fbc79e12b --- /dev/null +++ b/src/integrators/vi/position_momentum_midpoint.jl @@ -0,0 +1,28 @@ +@doc raw""" +Midpoint Variational Integrator in position-momentum form. + +We consider a discrete Lagrangian of the form +```math +L_d (q_{n}, q_{n+1}) = h \, L \bigg( \frac{q_{n} + q_{n+1}}{2}, \frac{q_{n+1} - q_{n}}{h} \bigg) , +``` +where $q_{n}$ approximates the solution $q(t_{n})$. +The Euler-Lagrange equations are computed as: +```math +\begin{aligned} +p_{n } &= - D_{1} L_{d} (q_{n}, q_{n+1}) = \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, \frac{q_{n+1} - q_{n}}{h} \bigg) , \\ +p_{n+1} &= \hphantom{-} D_{2} L_{d} (q_{n}, q_{n+1}) = \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, \frac{q_{n+1} - q_{n}}{h} \bigg) . +\end{aligned} +``` +The first equation can be solved implicitly for $q_{n+1}$ given $(q_{n}, p_{n})$. +The second equation can be used to explicitly compute $p_{n+1}$. +""" +struct PMVImidpoint <: VIMethod end + +isexplicit(method::PMVImidpoint) = false +isimplicit(method::PMVImidpoint) = true +issymmetric(method::PMVImidpoint) = missing +issymplectic(method::PMVImidpoint) = true + + +const PMVImidpointIntegrator{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:PMVImidpoint} + diff --git a/src/integrators/vi/position_momentum_trapezoidal.jl b/src/integrators/vi/position_momentum_trapezoidal.jl new file mode 100644 index 000000000..3dce35ec4 --- /dev/null +++ b/src/integrators/vi/position_momentum_trapezoidal.jl @@ -0,0 +1,28 @@ +@doc raw""" +Trapezoidal Variational Integrator in position-momentum form. + +We consider a discrete Lagrangian of the form +```math +L_d (q_{n}, q_{n+1}) = \frac{h}{2} \bigg[ L \bigg( q_{n}, \frac{q_{n+1} - q_{n}}{h} \bigg) + L \bigg( q_{n+1}, \frac{q_{n+1} - q_{n}}{h} \bigg) \bigg] , +``` +where $q_{n}$ approximates the solution $q(t_{n})$. +The Euler-Lagrange equations are computed as: +```math +\begin{aligned} +p_{n } &= - D_{1} L_{d} (q_{n}, q_{n+1}) = \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( q_{n}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{1}{2} \frac{\partial L}{\partial v} \bigg( q_{n}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{1}{2} \frac{\partial L}{\partial v} \bigg( q_{n+1}, \frac{q_{n+1} - q_{n}}{h} \bigg) , \\ +p_{n+1} &= \hphantom{-} D_{2} L_{d} (q_{n}, q_{n+1}) = \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( q_{n}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{1}{2} \frac{\partial L}{\partial v} \bigg( q_{n}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{1}{2} \frac{\partial L}{\partial v} \bigg( q_{n+1}, \frac{q_{n+1} - q_{n}}{h} \bigg) . +\end{aligned} +``` +The first equation can be solved implicitly for $q_{n+1}$ given $(q_{n}, p_{n})$. +The second equation can be used to explicitly compute $p_{n+1}$. +""" +struct PMVItrapezoidal <: VIMethod end + +isexplicit(method::PMVItrapezoidal) = false +isimplicit(method::PMVItrapezoidal) = true +issymmetric(method::PMVItrapezoidal) = missing +issymplectic(method::PMVItrapezoidal) = true + + +const PMVItrapezoidalIntegrator{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:PMVItrapezoidal} + diff --git a/src/integrators/vi/vi_methods.jl b/src/integrators/vi/vi_methods.jl new file mode 100644 index 000000000..628fbcfc0 --- /dev/null +++ b/src/integrators/vi/vi_methods.jl @@ -0,0 +1,5 @@ + +abstract type VIMethod <: LODEMethod end + +default_solver(::VIMethod) = Newton() +default_iguess(::VIMethod) = HermiteExtrapolation() diff --git a/src/integrators/vi/vprk_methods.jl b/src/integrators/vi/vprk_methods.jl index b7278835f..6d8b117b5 100644 --- a/src/integrators/vi/vprk_methods.jl +++ b/src/integrators/vi/vprk_methods.jl @@ -1,7 +1,7 @@ # Variational Partitioned Runge-Kutta Methods -abstract type VPRKMethod <: LODEMethod end +abstract type VPRKMethod <: VIMethod end GeometricBase.order(method::VPRKMethod) = RungeKutta.order(tableau(method)) diff --git a/test/integrators/hamilton_pontryagin_integrators_tests.jl b/test/integrators/hamilton_pontryagin_integrators_tests.jl new file mode 100644 index 000000000..e69de29bb diff --git a/test/runtests.jl b/test/runtests.jl index a74f05242..5e3e3d213 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,6 +12,8 @@ using SafeTestsets # @safetestset "Galerkin Variational Integrators " begin include("integrators/galerkin_integrators_tests.jl") end +@safetestset "Hamilton-Pontryagin Integrators " begin include("integrators/hamilton_pontryagin_integrators_tests.jl") end + @safetestset "Projection Methods " begin include("projections/projections_tests.jl") end @safetestset "Projection Methods with Implicit Equations " begin include("projections/projections_implicit_tests.jl") end @safetestset "Projection Methods with Variational Partitioned Runge-Kutta Integrators " begin include("projections/projections_vprk_tests.jl") end From b2692d717dd43c83fe7ee17a8d39cfce1adac4e6 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 29 Aug 2023 20:16:14 +0900 Subject: [PATCH 04/11] Implement position-momentum form of midpoint and trapezoidal variational integrator. --- src/Integrators.jl | 2 + src/integrators/vi/position_momentum_cache.jl | 41 ++++++++++++ .../vi/position_momentum_common.jl | 42 ++++++++++++ .../vi/position_momentum_midpoint.jl | 64 +++++++++++++++++- .../vi/position_momentum_trapezoidal.jl | 66 ++++++++++++++++++- src/integrators/vi/vi_methods.jl | 2 + src/integrators/vi/vprk_integrator.jl | 3 - src/integrators/vi/vprk_methods.jl | 2 - .../variational_integrators_tests.jl | 46 ++++++++++--- 9 files changed, 251 insertions(+), 17 deletions(-) create mode 100644 src/integrators/vi/position_momentum_cache.jl create mode 100644 src/integrators/vi/position_momentum_common.jl diff --git a/src/Integrators.jl b/src/Integrators.jl index 2ffac7bd2..96f7573e4 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -192,8 +192,10 @@ module Integrators export IntegratorVPRK include("integrators/vi/vi_methods.jl") + include("integrators/vi/position_momentum_cache.jl") include("integrators/vi/position_momentum_midpoint.jl") include("integrators/vi/position_momentum_trapezoidal.jl") + include("integrators/vi/position_momentum_common.jl") include("integrators/vi/vprk_methods.jl") include("integrators/vi/vprk_cache.jl") include("integrators/vi/vprk_integrator.jl") diff --git a/src/integrators/vi/position_momentum_cache.jl b/src/integrators/vi/position_momentum_cache.jl new file mode 100644 index 000000000..29df13af4 --- /dev/null +++ b/src/integrators/vi/position_momentum_cache.jl @@ -0,0 +1,41 @@ +@doc raw""" +Cache for variational integrator in position-momentum form. + +### Fields + +* `q`: internal stages of solution +* `Θ`: implicit function evaluated on solution +* `f`: vector field of implicit function +""" +struct IntegratorCachePMVI{DT,D} <: IODEIntegratorCache{DT,D} + x::Vector{DT} + + q::Vector{DT} + p::Vector{DT} + θ::Vector{DT} + f::Vector{DT} + + q̄::Vector{DT} + p̄::Vector{DT} + θ̄::Vector{DT} + f̄::Vector{DT} + + q̃::Vector{DT} + ṽ::Vector{DT} + θ̃::Vector{DT} + f̃::Vector{DT} + + function IntegratorCachePMVI{DT,D}() where {DT,D} + new(zeros(DT,D), + zeros(DT,D), zeros(DT,D), zeros(DT,D), zeros(DT,D), + zeros(DT,D), zeros(DT,D), zeros(DT,D), zeros(DT,D), + zeros(DT,D), zeros(DT,D), zeros(DT,D), zeros(DT,D)) + end +end + +function reset!(cache::IntegratorCachePMVI, t, q, p) + copyto!(cache.q̄, q) + copyto!(cache.p̄, p) +end + +nlsolution(cache::IntegratorCachePMVI) = cache.x diff --git a/src/integrators/vi/position_momentum_common.jl b/src/integrators/vi/position_momentum_common.jl new file mode 100644 index 000000000..a39953628 --- /dev/null +++ b/src/integrators/vi/position_momentum_common.jl @@ -0,0 +1,42 @@ + +function initial_guess!( + solstep::SolutionStepPODE{DT}, + problem::Union{IODEProblem,LODEProblem}, + method::Union{PMVImidpoint,PMVItrapezoidal}, + caches::CacheDict, + ::NonlinearSolver, + iguess::Union{InitialGuess,Extrapolation}) where {DT} + + # get cache and dimension + cache = caches[DT] + + # compute initial guess for solution q(n+1) + initialguess!(solstep.t, cache.q, cache.p, solstep, problem, iguess) + + # copy initial guess to solution vector + cache.x .= cache.q +end + +function integrate_step!( + solstep::SolutionStepPODE{DT,TT}, + problem::Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, + method::Union{PMVImidpoint,PMVItrapezoidal}, + caches::CacheDict, + solver::NonlinearSolver) where {DT,TT} + + # call nonlinear solver + solve!(caches[DT].x, (b,x) -> residual!(b, x, solstep, problem, method, caches), solver) + + # print solver status + # print_solver_status(int.solver.status, int.solver.params) + + # check if solution contains NaNs or error bounds are violated + # check_solver_status(int.solver.status, int.solver.params) + + # compute vector field at internal stages + components!(caches[DT].x, solstep, problem, method, caches) + + # compute final update + solstep.q .= caches[DT].q + solstep.p .= caches[DT].p +end diff --git a/src/integrators/vi/position_momentum_midpoint.jl b/src/integrators/vi/position_momentum_midpoint.jl index fbc79e12b..abf5f723b 100644 --- a/src/integrators/vi/position_momentum_midpoint.jl +++ b/src/integrators/vi/position_momentum_midpoint.jl @@ -20,9 +20,71 @@ struct PMVImidpoint <: VIMethod end isexplicit(method::PMVImidpoint) = false isimplicit(method::PMVImidpoint) = true -issymmetric(method::PMVImidpoint) = missing +issymmetric(method::PMVImidpoint) = true issymplectic(method::PMVImidpoint) = true const PMVImidpointIntegrator{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:PMVImidpoint} +function Cache{ST}(problem::Union{IODEProblem,LODEProblem}, method::PMVImidpoint; kwargs...) where {ST} + IntegratorCachePMVI{ST, ndims(problem)}(; kwargs...) +end + +@inline CacheType(ST, problem::Union{IODEProblem,LODEProblem}, ::PMVImidpoint) = IntegratorCachePMVI{ST, ndims(problem)} + +function Base.show(io::IO, int::PMVImidpointIntegrator) + print(io, "\nMidpoint variational integrator in position-momentum form with:\n") + print(io, " Timestep: $(timestep(int))\n") +end + + +function components!( + x::Vector{ST}, + solstep::SolutionStepPODE{DT,TT}, + problem::Union{IODEProblem,LODEProblem}, + method::PMVImidpoint, + caches::CacheDict) where {ST,DT,TT} + + # get cache and dimension + cache = caches[ST] + D = ndims(problem) + + # set some local variables for convenience and clarity + local t̃ = solstep.t̄ + timestep(problem) / 2 + + # copy x to q + cache.q .= x[1:D] + + # compute q̃ and ṽ + cache.q̃ .= (cache.q .+ cache.q̄) ./ 2 + cache.ṽ .= (cache.q .- cache.q̄) ./ timestep(problem) + + # compute Θ̃ = ϑ(q̃,ṽ) and f̃ = f(q̃,ṽ) + functions(problem).ϑ(cache.θ̃, t̃, cache.q̃, cache.ṽ) + functions(problem).f(cache.f̃, t̃, cache.q̃, cache.ṽ) + + # compute p + cache.p .= cache.θ̃ .+ timestep(problem) .* cache.f̃ ./ 2 +end + + +function residual!( + b::Vector{ST}, + x::Vector{ST}, + solstep::SolutionStepPODE, + problem::Union{IODEProblem,LODEProblem}, + method::PMVImidpoint, + caches::CacheDict) where {ST} + + # get cache for internal stages + cache = caches[ST] + + # copy previous solution from solstep to cache + reset!(cache, current(solstep)...) + + # compute stages from nonlinear solver solution x + components!(x, solstep, problem, method, caches) + + # compute b + b .= cache.θ̃ .- cache.p̄ .- timestep(problem) .* cache.f̃ ./ 2 +end diff --git a/src/integrators/vi/position_momentum_trapezoidal.jl b/src/integrators/vi/position_momentum_trapezoidal.jl index 3dce35ec4..b93d86bad 100644 --- a/src/integrators/vi/position_momentum_trapezoidal.jl +++ b/src/integrators/vi/position_momentum_trapezoidal.jl @@ -20,9 +20,73 @@ struct PMVItrapezoidal <: VIMethod end isexplicit(method::PMVItrapezoidal) = false isimplicit(method::PMVItrapezoidal) = true -issymmetric(method::PMVItrapezoidal) = missing +issymmetric(method::PMVItrapezoidal) = true issymplectic(method::PMVItrapezoidal) = true const PMVItrapezoidalIntegrator{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:PMVItrapezoidal} +function Cache{ST}(problem::Union{IODEProblem,LODEProblem}, method::PMVItrapezoidal; kwargs...) where {ST} + IntegratorCachePMVI{ST, ndims(problem)}(; kwargs...) +end + +@inline CacheType(ST, problem::Union{IODEProblem,LODEProblem}, ::PMVItrapezoidal) = IntegratorCachePMVI{ST, ndims(problem)} + +function Base.show(io::IO, int::PMVItrapezoidalIntegrator) + print(io, "\nTrapezoidal variational integrator in position-momentum form with:\n") + print(io, " Timestep: $(timestep(int))\n") +end + + +function components!( + x::Vector{ST}, + solstep::SolutionStepPODE{DT,TT}, + problem::Union{IODEProblem,LODEProblem}, + method::PMVItrapezoidal, + caches::CacheDict) where {ST,DT,TT} + + # get cache and dimension + cache = caches[ST] + D = ndims(problem) + + # set some local variables for convenience and clarity + local t̄ = solstep.t̄ + local t = solstep.t̄ + timestep(problem) + + # copy x to q + cache.q .= x[1:D] + + # compute v + cache.ṽ .= (cache.q .- cache.q̄) ./ timestep(problem) + + # compute Θ = ϑ(q,ṽ) and f = f(q,ṽ) + functions(problem).ϑ(cache.θ̄, t̄, cache.q̄, cache.ṽ) + functions(problem).ϑ(cache.θ, t, cache.q, cache.ṽ) + functions(problem).f(cache.f̄, t̄, cache.q̄, cache.ṽ) + functions(problem).f(cache.f, t, cache.q, cache.ṽ) + + # compute p + cache.p .= (cache.θ .+ cache.θ̄) ./ 2 .+ timestep(problem) .* cache.f ./ 2 +end + + +function residual!( + b::Vector{ST}, + x::Vector{ST}, + solstep::SolutionStepPODE, + problem::Union{IODEProblem,LODEProblem}, + method::PMVItrapezoidal, + caches::CacheDict) where {ST} + + # get cache for internal stages + cache = caches[ST] + + # copy previous solution from solstep to cache + reset!(cache, current(solstep)...) + + # compute stages from nonlinear solver solution x + components!(x, solstep, problem, method, caches) + + # compute b + b .= (cache.θ .+ cache.θ̄) ./ 2 .- cache.p̄ .- timestep(problem) .* cache.f̄ ./ 2 +end diff --git a/src/integrators/vi/vi_methods.jl b/src/integrators/vi/vi_methods.jl index 628fbcfc0..379d3c8c6 100644 --- a/src/integrators/vi/vi_methods.jl +++ b/src/integrators/vi/vi_methods.jl @@ -1,5 +1,7 @@ abstract type VIMethod <: LODEMethod end +isiodemethod(::Union{VIMethod, Type{<:VIMethod}}) = true + default_solver(::VIMethod) = Newton() default_iguess(::VIMethod) = HermiteExtrapolation() diff --git a/src/integrators/vi/vprk_integrator.jl b/src/integrators/vi/vprk_integrator.jl index 2bc502e41..553a15a11 100644 --- a/src/integrators/vi/vprk_integrator.jl +++ b/src/integrators/vi/vprk_integrator.jl @@ -25,9 +25,6 @@ const IntegratorVPRK{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LOD description(::IntegratorVPRK) = "Variational Partitioned Runge-Kutta Integrator" -default_solver(::VPRKMethod) = Newton() -default_iguess(::VPRKMethod) = HermiteExtrapolation() - initmethod(method::VPRKMethod) = VPRK(method) solversize(problem::VPRKProblem, method::VPRKMethod) = diff --git a/src/integrators/vi/vprk_methods.jl b/src/integrators/vi/vprk_methods.jl index 6d8b117b5..8eba13a42 100644 --- a/src/integrators/vi/vprk_methods.jl +++ b/src/integrators/vi/vprk_methods.jl @@ -8,8 +8,6 @@ GeometricBase.order(method::VPRKMethod) = RungeKutta.order(tableau(method)) @inline nstages(method::VPRKMethod) = nstages(tableau(method)) @inline eachstage(method::VPRKMethod) = eachstage(tableau(method)) -isiodemethod(::Union{VPRKMethod, Type{<:VPRKMethod}}) = true - isexplicit(method::VPRKMethod) = RungeKutta.isexplicit(tableau(method)) isimplicit(method::VPRKMethod) = RungeKutta.isimplicit(tableau(method)) issymmetric(method::VPRKMethod) = RungeKutta.issymmetric(tableau(method)) diff --git a/test/integrators/variational_integrators_tests.jl b/test/integrators/variational_integrators_tests.jl index 7074a560b..50ccade85 100644 --- a/test/integrators/variational_integrators_tests.jl +++ b/test/integrators/variational_integrators_tests.jl @@ -3,69 +3,95 @@ using GeometricProblems.HarmonicOscillator using Test -iode = iodeproblem() +lode = lodeproblem() pref = exact_solution(podeproblem()) +@testset "$(rpad("Vartiational integrators",80))" begin + + sol = integrate(lode, PMVImidpoint()) + @test relative_maximum_error(sol.q, pref.q) < 4E-4 + + sol = integrate(lode, PMVItrapezoidal()) + @test relative_maximum_error(sol.q, pref.q) < 4E-4 + +end + + @testset "$(rpad("Vartiational Partitioned Runge-Kutta integrators",80))" begin - sol = integrate(iode, VPRKGauss(1)) + sol = integrate(lode, VPRKGauss(1)) # println(relative_maximum_error(sol.q, pref.q)) # println(relative_maximum_error(sol.p, pref.p)) # println() @test relative_maximum_error(sol.q, pref.q) < 4E-4 # @test relative_maximum_error(sol.p, pref.p) < 4E-4 - sol = integrate(iode, VPRKGauss(2)) + ref = integrate(lode, PMVImidpoint()) + @test relative_maximum_error(sol.q, ref.q) < 8*eps() + + sol = integrate(lode, VPRKGauss(2)) # println(relative_maximum_error(sol.q, pref.q)) # println(relative_maximum_error(sol.p, pref.p)) # println() @test relative_maximum_error(sol.q, pref.q) < 4E-8 # @test relative_maximum_error(sol.p, pref.p) < 4E-8 - sol = integrate(iode, VPRKGauss(3)) + sol = integrate(lode, VPRKGauss(3)) # println(relative_maximum_error(sol.q, pref.q)) # println(relative_maximum_error(sol.p, pref.p)) # println() @test relative_maximum_error(sol.q, pref.q) < 8E-13 # @test relative_maximum_error(sol.p, pref.p) < 8E-13 - sol = integrate(iode, VPRKLobattoIIIAIIIĀ(2)) + sol = integrate(lode, VPRKGauss(4)) + # println(relative_maximum_error(sol.q, pref.q)) + # println(relative_maximum_error(sol.p, pref.p)) + # println() + @test relative_maximum_error(sol.q, pref.q) < 2E-16 + # @test relative_maximum_error(sol.p, pref.p) < 2E-16 + + + sol = integrate(lode, VPRKLobattoIIIAIIIĀ(2)) # println(relative_maximum_error(sol.q, pref.q)) # println(relative_maximum_error(sol.p, pref.p)) # println() @test relative_maximum_error(sol.q, pref.q) < 2E-4 # @test relative_maximum_error(sol.p, pref.p) < 2E-4 - sol = integrate(iode, VPRKLobattoIIIAIIIĀ(3)) + ref = integrate(lode, PMVItrapezoidal()) + @test relative_maximum_error(sol.q, ref.q) < 8*eps() + + sol = integrate(lode, VPRKLobattoIIIAIIIĀ(3)) # println(relative_maximum_error(sol.q, pref.q)) # println(relative_maximum_error(sol.p, pref.p)) # println() @test relative_maximum_error(sol.q, pref.q) < 8E-9 # @test relative_maximum_error(sol.p, pref.p) < 8E-9 - sol = integrate(iode, VPRKLobattoIIIAIIIĀ(4)) + sol = integrate(lode, VPRKLobattoIIIAIIIĀ(4)) # println(relative_maximum_error(sol.q, pref.q)) # println(relative_maximum_error(sol.p, pref.p)) # println() @test relative_maximum_error(sol.q, pref.q) < 2E-13 # @test relative_maximum_error(sol.p, pref.p) < 2E-13 - sol = integrate(iode, VPRKLobattoIIIBIIIB̄(2)) + + sol = integrate(lode, VPRKLobattoIIIBIIIB̄(2)) # println(relative_maximum_error(sol.q, pref.q)) # println(relative_maximum_error(sol.p, pref.p)) # println() @test relative_maximum_error(sol.q, pref.q) < 4E-4 # @test relative_maximum_error(sol.p, pref.p) < 4E-4 - sol = integrate(iode, VPRKLobattoIIIBIIIB̄(3)) + sol = integrate(lode, VPRKLobattoIIIBIIIB̄(3)) # println(relative_maximum_error(sol.q, pref.q)) # println(relative_maximum_error(sol.p, pref.p)) # println() @test relative_maximum_error(sol.q, pref.q) < 4E-8 # @test relative_maximum_error(sol.p, pref.p) < 4E-8 - sol = integrate(iode, VPRKLobattoIIIBIIIB̄(4)) + sol = integrate(lode, VPRKLobattoIIIBIIIB̄(4)) # println(relative_maximum_error(sol.q, pref.q)) # println(relative_maximum_error(sol.p, pref.p)) # println() From f8fe5d1494fe04874b6b2ae966bb68534b454496 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 29 Aug 2023 20:41:14 +0900 Subject: [PATCH 05/11] Minor fixes and cleanup in DVIs. --- src/integrators/dvi/integrators_cmdvi.jl | 6 +----- src/integrators/dvi/integrators_ctdvi.jl | 6 +----- src/integrators/dvi/integrators_dvi_a.jl | 6 +----- src/integrators/dvi/integrators_dvi_b.jl | 6 +----- src/integrators/dvi/integrators_dvrk.jl | 6 +----- src/integrators/dvi/methods.jl | 3 +++ 6 files changed, 8 insertions(+), 25 deletions(-) diff --git a/src/integrators/dvi/integrators_cmdvi.jl b/src/integrators/dvi/integrators_cmdvi.jl index 85d4f1d79..ad91dbd22 100644 --- a/src/integrators/dvi/integrators_cmdvi.jl +++ b/src/integrators/dvi/integrators_cmdvi.jl @@ -8,7 +8,7 @@ Degenerate variational integrator cache. * `Θ`: implicit function evaluated on solution * `f`: vector field of implicit function """ -struct IntegratorCacheCMDVI{DT,D} <: ODEIntegratorCache{DT,D} +struct IntegratorCacheCMDVI{DT,D} <: IODEIntegratorCache{DT,D} x::Vector{DT} q::Vector{DT} @@ -42,10 +42,6 @@ Midpoint Degenerate Variational Integrator. const IntegratorCMDVI{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:CMDVI} -default_solver(::CMDVI) = Newton() -default_iguess(::CMDVI) = HermiteExtrapolation() - - function Base.show(io::IO, int::IntegratorCMDVI) print(io, "\nMidpoint Degenerate Variational Integrator with:\n") print(io, " Timestep: $(timestep(int))\n") diff --git a/src/integrators/dvi/integrators_ctdvi.jl b/src/integrators/dvi/integrators_ctdvi.jl index a1d66f492..f0f3cfbe0 100644 --- a/src/integrators/dvi/integrators_ctdvi.jl +++ b/src/integrators/dvi/integrators_ctdvi.jl @@ -8,7 +8,7 @@ Degenerate variational integrator cache. * `Θ`: implicit function evaluated on solution * `f`: vector field of implicit function """ -struct IntegratorCacheCTDVI{DT,D} <: ODEIntegratorCache{DT,D} +struct IntegratorCacheCTDVI{DT,D} <: IODEIntegratorCache{DT,D} x::Vector{DT} q::Vector{DT} @@ -47,10 +47,6 @@ Trapezoidal Degenerate Variational Integrator. const IntegratorCTDVI{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:CTDVI} -default_solver(::CTDVI) = Newton() -default_iguess(::CTDVI) = HermiteExtrapolation() - - function Base.show(io::IO, int::IntegratorCTDVI) print(io, "\nTrapezoidal Degenerate Variational Integrator with:\n") print(io, " Timestep: $(timestep(int))\n") diff --git a/src/integrators/dvi/integrators_dvi_a.jl b/src/integrators/dvi/integrators_dvi_a.jl index 0783cb9c5..56eea058a 100644 --- a/src/integrators/dvi/integrators_dvi_a.jl +++ b/src/integrators/dvi/integrators_dvi_a.jl @@ -8,7 +8,7 @@ Degenerate variational integrator cache. * `Θ`: implicit function evaluated on solution * `f`: vector field of implicit function """ -struct IntegratorCacheDVIA{DT,D} <: ODEIntegratorCache{DT,D} +struct IntegratorCacheDVIA{DT,D} <: IODEIntegratorCache{DT,D} x::Vector{DT} q::Vector{DT} @@ -42,10 +42,6 @@ Symplectic Euler-A Degenerate Variational Integrator. const IntegratorDVIA{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:DVIA} -default_solver(::DVIA) = Newton() -default_iguess(::DVIA) = HermiteExtrapolation() - - function Base.show(io::IO, int::IntegratorDVIA) print(io, "\nDegenerate Variational Integrator (Euler-A) with:\n") print(io, " Timestep: $(timestep(int))\n") diff --git a/src/integrators/dvi/integrators_dvi_b.jl b/src/integrators/dvi/integrators_dvi_b.jl index ae156264a..c05c75478 100644 --- a/src/integrators/dvi/integrators_dvi_b.jl +++ b/src/integrators/dvi/integrators_dvi_b.jl @@ -8,7 +8,7 @@ Degenerate variational integrator cache. * `Θ`: implicit function evaluated on solution * `f`: vector field of implicit function """ -struct IntegratorCacheDVIB{DT,D} <: ODEIntegratorCache{DT,D} +struct IntegratorCacheDVIB{DT,D} <: IODEIntegratorCache{DT,D} x::Vector{DT} q::Vector{DT} @@ -42,10 +42,6 @@ Symplectic Euler-B Degenerate Variational Integrator. const IntegratorDVIB{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:DVIB} -default_solver(::DVIB) = Newton() -default_iguess(::DVIB) = HermiteExtrapolation() - - function Base.show(io::IO, int::IntegratorDVIB) print(io, "\nDegenerate Variational Integrator (Euler-B) with:\n") print(io, " Timestep: $(timestep(int))\n") diff --git a/src/integrators/dvi/integrators_dvrk.jl b/src/integrators/dvi/integrators_dvrk.jl index a9dab6640..fd286ab31 100644 --- a/src/integrators/dvi/integrators_dvrk.jl +++ b/src/integrators/dvi/integrators_dvrk.jl @@ -65,7 +65,7 @@ Degenerate Variational Runge-Kutta integrator cache. * `Θ`: implicit function of internal stages * `F`: vector field of implicit function """ -struct IntegratorCacheDVRK{DT,D,S} <: ODEIntegratorCache{DT,D} +struct IntegratorCacheDVRK{DT,D,S} <: IODEIntegratorCache{DT,D} x::Vector{DT} q::Vector{DT} @@ -99,10 +99,6 @@ end const IntegratorDVRK{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:DVRK} -default_solver(::DVRK) = Newton() -default_iguess(::DVRK) = HermiteExtrapolation() - - function Base.show(io::IO, int::IntegratorDVRK) print(io, "\nRunge-Kutta Integrator for Degenerate Lagrangians with:\n") print(io, " Timestep: $(timestep(int))\n") diff --git a/src/integrators/dvi/methods.jl b/src/integrators/dvi/methods.jl index bca6b5d47..8bf3f3209 100644 --- a/src/integrators/dvi/methods.jl +++ b/src/integrators/dvi/methods.jl @@ -5,6 +5,9 @@ issymplectic(::Union{DVIMethod, Type{<:DVIMethod}}) = true isexplicit(::Union{DVIMethod, Type{<:DVIMethod}}) = false isimplicit(::Union{DVIMethod, Type{<:DVIMethod}}) = true +default_solver(::DVIMethod) = Newton() +default_iguess(::DVIMethod) = HermiteExtrapolation() + "Symplectic Euler-A Degenerate Variational Integrator." struct DVIA <: DVIMethod end From d48b40a80019b26b59dada353f421182f4c21411 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Wed, 30 Aug 2023 00:27:30 +0900 Subject: [PATCH 06/11] Implement trapezoidal Hamilton-Pontryagin integrator. --- src/Integrators.jl | 2 + src/integrators/hpi/hpi_methods.jl | 5 + src/integrators/hpi/hpi_midpoint.jl | 8 +- src/integrators/hpi/hpi_trapezoidal.jl | 99 ++++++++++++++++++- .../hamilton_pontryagin_integrators_tests.jl | 46 +++++++++ 5 files changed, 154 insertions(+), 6 deletions(-) diff --git a/src/Integrators.jl b/src/Integrators.jl index 96f7573e4..2c8119b7a 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -227,8 +227,10 @@ module Integrators export HPImidpointIntegrator, HPItrapezoidalIntegrator include("integrators/hpi/hpi_methods.jl") + include("integrators/hpi/hpi_cache.jl") include("integrators/hpi/hpi_midpoint.jl") include("integrators/hpi/hpi_trapezoidal.jl") + include("integrators/hpi/hpi_common.jl") export NoProjection, projection diff --git a/src/integrators/hpi/hpi_methods.jl b/src/integrators/hpi/hpi_methods.jl index bd717e6fd..ad63370ad 100644 --- a/src/integrators/hpi/hpi_methods.jl +++ b/src/integrators/hpi/hpi_methods.jl @@ -1,2 +1,7 @@ abstract type HPIMethod <: LODEMethod end + +isiodemethod(::Union{HPIMethod, Type{<:HPIMethod}}) = true + +default_solver(::HPIMethod) = Newton() +default_iguess(::HPIMethod) = HermiteExtrapolation() diff --git a/src/integrators/hpi/hpi_midpoint.jl b/src/integrators/hpi/hpi_midpoint.jl index a6da90c7c..8195df2d2 100644 --- a/src/integrators/hpi/hpi_midpoint.jl +++ b/src/integrators/hpi/hpi_midpoint.jl @@ -8,7 +8,7 @@ L_d (q_{n}, q_{n+1}) = h \, L \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) where $q_{n}$ approximates the solution $q(t_n)$ and $v_{n+1/2}$ is the velocity, which is assumed to be constant in the interval $[t_{n}, t_{n+1}]$. The discrete Hamilton-Pontryagin action reads ```math -A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_{d}^{\alpha,\beta} (q_{n}, q_{n+1}, v_{n+1/2}) + \left< p_{n+1/2} , v_{n+1/2} - \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \right> \bigg] , +A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_{d}^{\alpha,\beta} (q_{n}, q_{n+1}, v_{n+1/2}) + \left< p_{n+1/2} , \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) - v_{n+1/2} \right> \bigg] , ``` where $\phi_h$ is a map that computes the velocity $v_{n+1/2}$ as a function of $q_{n}$, $q_{n+1}$ and a set of parameters $a_{n,n+1}$. A trivial example of such a map that does not depend on any parameters $a_{n,n+1}$ is @@ -22,8 +22,8 @@ The equations of motion, that are solved by this integrator, is computed as: \begin{aligned} 0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n-1} + q_{n}}{2}, v_{n-1/2} \bigg) + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) \\ - &- h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} - - h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} , \\ + &+ h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + + h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} , \\ 0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} , \\ p_{n+1/2} &= \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) , \\ v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) . @@ -41,6 +41,8 @@ struct HPImidpoint{ϕT, D₁ϕT, D₂ϕT, DₐϕT, PT} <: HPIMethod params::PT end +nparams(method::HPImidpoint) = length(method.params) + isexplicit(method::HPImidpoint) = false isimplicit(method::HPImidpoint) = true issymmetric(method::HPImidpoint) = missing diff --git a/src/integrators/hpi/hpi_trapezoidal.jl b/src/integrators/hpi/hpi_trapezoidal.jl index 63dbe226d..932852dde 100644 --- a/src/integrators/hpi/hpi_trapezoidal.jl +++ b/src/integrators/hpi/hpi_trapezoidal.jl @@ -8,7 +8,7 @@ L_d (q_{n}, q_{n+1}) = \frac{h}{2} \big[ L (q_{n}, v_{n+1/2}) + L (q_{n+1}, v_{n where $q_{n}$ approximates the solution $q(t_n)$ and $v_{n+1/2}$ is the velocity, which is assumed to be constant in the interval $[t_{n}, t_{n+1}]$. The discrete Hamilton-Pontryagin action reads ```math -A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_d (q_{n}, q_{n+1}) + h \left< p_{n+1} , v_{n+1/2} - \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \right> \bigg] , +A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_d (q_{n}, q_{n+1}) + h \left< p_{n+1} , \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) - v_{n+1/2} \right> \bigg] , ``` where $\phi_h$ is a map that computes the velocity $v_{n+1/2}$ as a function of $q_{n}$, $q_{n+1}$ and a set of parameters $a_{n,n+1}$. A trivial example of such a map that does not depend on any parameters $a_{n,n+1}$ is @@ -22,8 +22,8 @@ The equations of motion, that are solved by this integrator, is computed as: \begin{aligned} 0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2}) + \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n+1}, v_{n+1/2}) \\ - &- h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} - - h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n} , \\ + &+ h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} + + h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n} , \\ 0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} , \\ p_{n+1} &= \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n}, v_{n+1/2}) + \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n+1}, v_{n+1/2}) , \\ @@ -42,6 +42,8 @@ struct HPItrapezoidal{ϕT, D₁ϕT, D₂ϕT, DₐϕT, PT} <: HPIMethod params::PT end +nparams(method::HPItrapezoidal) = length(method.params) + isexplicit(method::HPItrapezoidal) = false isimplicit(method::HPItrapezoidal) = true issymmetric(method::HPItrapezoidal) = missing @@ -49,3 +51,94 @@ issymplectic(method::HPItrapezoidal) = true const HPItrapezoidalIntegrator{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:HPItrapezoidal} + +function Cache{ST}(problem::Union{IODEProblem,LODEProblem}, method::HPItrapezoidal; kwargs...) where {ST} + IntegratorCacheHPI{ST, ndims(problem), nparams(method)}(; kwargs...) +end + +@inline CacheType(ST, problem::Union{IODEProblem,LODEProblem}, method::HPItrapezoidal) = IntegratorCacheHPI{ST, ndims(problem), nparams(method)} + +function Base.show(io::IO, int::HPItrapezoidalIntegrator) + print(io, "\nHamilton-Pontryagin Integrator using trapezoidal quadrature with:\n") + print(io, " Timestep: $(timestep(int))\n") +end + + +function components!( + x::Vector{ST}, + solstep::SolutionStepPODE{DT,TT}, + problem::Union{IODEProblem,LODEProblem}, + method::HPItrapezoidal, + caches::CacheDict) where {ST,DT,TT} + + # get cache and dimension + cache = caches[ST] + D = ndims(problem) + A = nparams(method) + + # set some local variables for convenience and clarity + local t̄ = solstep.t̄ + local t = solstep.t̄ + timestep(problem) + + # copy x to q + cache.q .= x[1:D] + cache.a .= x[D+1:D+A] + + # compute v + method.ϕ(cache.ṽ, cache.q̄, cache.q, cache.a, timestep(problem)) + + # compute Θ = ϑ(q,ṽ) and f = f(q,ṽ) + functions(problem).ϑ(cache.θ̄, t̄, cache.q̄, cache.ṽ) + functions(problem).ϑ(cache.θ, t, cache.q, cache.ṽ) + functions(problem).f(cache.f̄, t̄, cache.q̄, cache.ṽ) + functions(problem).f(cache.f, t, cache.q, cache.ṽ) + + # compute derivatives of ϕ + method.D₁ϕ(cache.D₁ϕ, cache.q̄, cache.q, cache.a, timestep(problem)) + method.D₂ϕ(cache.D₂ϕ, cache.q̄, cache.q, cache.a, timestep(problem)) + method.Dₐϕ(cache.Dₐϕ, cache.q̄, cache.q, cache.a, timestep(problem)) + + # compute p + cache.θ̃ .= (cache.θ .+ cache.θ̄) ./ 2 + cache.p .= timestep(problem) .* cache.f ./ 2 + for i in 1:D + for j in 1:D + cache.p[i] += timestep(problem) * cache.D₂ϕ[i,j] * cache.θ̃[j] + end + end +end + + +function residual!( + b::Vector{ST}, + x::Vector{ST}, + solstep::SolutionStepPODE, + problem::Union{IODEProblem,LODEProblem}, + method::HPItrapezoidal, + caches::CacheDict) where {ST} + + # get cache for internal stages + cache = caches[ST] + D = ndims(problem) + A = nparams(method) + + # copy previous solution from solstep to cache + reset!(cache, current(solstep)...) + + # compute stages from nonlinear solver solution x + components!(x, solstep, problem, method, caches) + + # compute b + for i in 1:D + b[i] = cache.p̄[i] + timestep(problem) * cache.f̄[i] / 2 + for j in 1:D + b[i] += timestep(problem) * cache.D₁ϕ[i,j] * cache.θ̃[j] + end + end + for i in 1:A + b[D+i] = 0 + for j in 1:D + b[D+i] += cache.Dₐϕ[i,j] * cache.θ̃[j] + end + end +end diff --git a/test/integrators/hamilton_pontryagin_integrators_tests.jl b/test/integrators/hamilton_pontryagin_integrators_tests.jl index e69de29bb..17d0c61d3 100644 --- a/test/integrators/hamilton_pontryagin_integrators_tests.jl +++ b/test/integrators/hamilton_pontryagin_integrators_tests.jl @@ -0,0 +1,46 @@ +using GeometricIntegrators +using GeometricProblems.HarmonicOscillator +using Test + + +lode = lodeproblem() +pref = exact_solution(podeproblem()) + + +@testset "$(rpad("Hamilton-Pontryagin integrators",80))" begin + + ϕ(v, q̄, q, a, Δt) = v .= (q .- q̄) ./ Δt + + function D₁ϕ(d, q̄, q, a, Δt) + d .= 0 + for i in eachindex(q) + d[i,i] = - 1 / Δt + end + end + + function D₂ϕ(d, q̄, q, a, Δt) + d .= 0 + for i in eachindex(q) + d[i,i] = + 1 / Δt + end + end + + function Dₐϕ(d, q̄, q, a, Δt) + d .= 0 + end + + + # sol = integrate(lode, HPImidpoint(ϕ, D₁ϕ, D₂ϕ, Dₐϕ, Float64[])) + # @test relative_maximum_error(sol.q, pref.q) < 4E-4 + + # ref = integrate(lode, PMVImidpoint()) + # @test relative_maximum_error(sol.q, ref.q) < 8*eps() + + + sol = integrate(lode, HPItrapezoidal(ϕ, D₁ϕ, D₂ϕ, Dₐϕ, Float64[])) + @test relative_maximum_error(sol.q, pref.q) < 4E-4 + + ref = integrate(lode, PMVItrapezoidal()) + @test relative_maximum_error(sol.q, ref.q) < 8*eps() + +end From ceae76ce76aff372fb17f2355c5cb5821fbf2e3d Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Wed, 30 Aug 2023 00:27:47 +0900 Subject: [PATCH 07/11] Minor refactoring in position-momentum variational integrators. --- src/integrators/vi/position_momentum_midpoint.jl | 2 +- src/integrators/vi/position_momentum_trapezoidal.jl | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/integrators/vi/position_momentum_midpoint.jl b/src/integrators/vi/position_momentum_midpoint.jl index abf5f723b..a325a13fd 100644 --- a/src/integrators/vi/position_momentum_midpoint.jl +++ b/src/integrators/vi/position_momentum_midpoint.jl @@ -64,7 +64,7 @@ function components!( functions(problem).f(cache.f̃, t̃, cache.q̃, cache.ṽ) # compute p - cache.p .= cache.θ̃ .+ timestep(problem) .* cache.f̃ ./ 2 + cache.p .= cache.p̄ .+ timestep(problem) .* cache.f̃ end diff --git a/src/integrators/vi/position_momentum_trapezoidal.jl b/src/integrators/vi/position_momentum_trapezoidal.jl index b93d86bad..adfd14582 100644 --- a/src/integrators/vi/position_momentum_trapezoidal.jl +++ b/src/integrators/vi/position_momentum_trapezoidal.jl @@ -66,7 +66,8 @@ function components!( functions(problem).f(cache.f, t, cache.q, cache.ṽ) # compute p - cache.p .= (cache.θ .+ cache.θ̄) ./ 2 .+ timestep(problem) .* cache.f ./ 2 + cache.θ̃ .= (cache.θ .+ cache.θ̄) ./ 2 + cache.p .= cache.p̄ .+ timestep(problem) .* (cache.f .+ cache.f̄) ./ 2 end @@ -88,5 +89,5 @@ function residual!( components!(x, solstep, problem, method, caches) # compute b - b .= (cache.θ .+ cache.θ̄) ./ 2 .- cache.p̄ .- timestep(problem) .* cache.f̄ ./ 2 + b .= cache.θ̃ .- cache.p̄ .- timestep(problem) .* cache.f̄ ./ 2 end From 6bfbf6436a20d4a2cd1cd620683d100f0ed899a9 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Wed, 30 Aug 2023 11:10:01 +0900 Subject: [PATCH 08/11] Update documentation for trapezoidal Hamilton-Pontryagin integrator. --- docs/make.jl | 4 +- docs/src/integrators/hpi.md | 66 ++++++++++++++++++++++++++ src/integrators/hpi/hpi_trapezoidal.jl | 27 ++++++----- 3 files changed, 83 insertions(+), 14 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 3a51e7fe3..d091427be 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -55,8 +55,8 @@ makedocs(bib, "DVI" => "integrators/dvi.md", # "CGVI" => "integrators/cgvi.md", # "DGVI" => "integrators/dgvi.md", - # "HPI" => "integrators/hpi.md", - # "HPG" => "integrators/hpg.md", + "Hamilton-Pontryagin" => "integrators/hpi.md", + # "Hamilton-Pontryagin-Galerkin" => "integrators/hpg.md", ], "Modules" => [ # "Discontinuities" => "modules/discontinuities.md", diff --git a/docs/src/integrators/hpi.md b/docs/src/integrators/hpi.md index 7a2b6074b..ae98b6278 100644 --- a/docs/src/integrators/hpi.md +++ b/docs/src/integrators/hpi.md @@ -1 +1,67 @@ # [Hamilton-Pontryagin Integrators](@id hpi) + + + +## Trapezoidal Discrete Lagrangian + +We consider a discrete Lagrangian of the form +```math +L_d (q_{n}, q_{n+1}) = \frac{h}{2} \big[ L (q_{n}, v_{n+1/2}) + L (q_{n+1}, v_{n+1/2}) \big] , +``` +where $q_{n}$ approximates the solution $q(t_n)$ and $v_{n+1/2}$ is the velocity, which is assumed to be constant in the interval $[t_{n}, t_{n+1}]$. +The discrete Hamilton-Pontryagin action reads +```math +A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_d (q_{n}, q_{n+1}) + h \left< p_{n+1/2} , \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) - v_{n+1/2} \right> \bigg] , +``` +where $\phi_h$ is a map that computes the velocity $v_{n+1/2}$ as a function of $q_{n}$, $q_{n+1}$ and a set of parameters $a_{n,n+1}$. +A trivial example of such a map that does not depend on any parameters $a_{n,n+1}$ is +```math +\phi_h (q_{n}, q_{n+1}; a_{n,n+1}) = \frac{q_{n+1} - q_{n}}{h} . +``` +In order to solve the discrete Euler-Lagrange equations, the user needs to specify the map $\phi_h$, its gradients with respect to $q_{n}$ and $q_{n+1}$, denoted by $D_1 \phi_h$ and $D_2 \phi_h$, respectively, the gradient with respect to the parameters, denoted by $D_a \phi_h$, and an initial set of parameters $a_{0}$. + +The equations of motion, that are solved by this integrator, are computed as: +```math +\begin{aligned} +0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2}) + + \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n-1/2}) \\ + &+ h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + + h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} , \\ +0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ +v_{n+1/2} +&= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\ +p_{n+1/2} +&= \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n}, v_{n+1/2}) + + \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n+1}, v_{n+1/2}) . +\end{aligned} +``` + +Upon defining the momentum +```math +p_{n} += h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} ++ \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n-1/2}) +``` +we can rewrite the equations of motion as +```math +\begin{aligned} +0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2}) + + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + + p_{n} , \\ +0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ +v_{n+1/2} +&= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\ +p_{n+1/2} +&= \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n}, v_{n+1/2}) + + \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n+1}, v_{n+1/2}) , \\ +p_{n+1} +&= h \, D_2 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + + \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n+1}, v_{n+1/2}) . +\end{aligned} +``` +Given $(q_{n}, p_{n})$, the first four equations can be solved for $q_{n+1}$, where $v_{n+1/2}$, $p_{n+1/2}$, and $a_{n,n+1}$ are treated as internal variables similar to the internal stages of a Runge-Kutta method, and the last equation provides an explicit update for $p_{n+1}$. + + +## Midpoint Discrete Lagrangian + + diff --git a/src/integrators/hpi/hpi_trapezoidal.jl b/src/integrators/hpi/hpi_trapezoidal.jl index 932852dde..1ebe772b0 100644 --- a/src/integrators/hpi/hpi_trapezoidal.jl +++ b/src/integrators/hpi/hpi_trapezoidal.jl @@ -8,7 +8,7 @@ L_d (q_{n}, q_{n+1}) = \frac{h}{2} \big[ L (q_{n}, v_{n+1/2}) + L (q_{n+1}, v_{n where $q_{n}$ approximates the solution $q(t_n)$ and $v_{n+1/2}$ is the velocity, which is assumed to be constant in the interval $[t_{n}, t_{n+1}]$. The discrete Hamilton-Pontryagin action reads ```math -A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_d (q_{n}, q_{n+1}) + h \left< p_{n+1} , \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) - v_{n+1/2} \right> \bigg] , +A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_d (q_{n}, q_{n+1}) + h \left< p_{n+1/2} , \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) - v_{n+1/2} \right> \bigg] , ``` where $\phi_h$ is a map that computes the velocity $v_{n+1/2}$ as a function of $q_{n}$, $q_{n+1}$ and a set of parameters $a_{n,n+1}$. A trivial example of such a map that does not depend on any parameters $a_{n,n+1}$ is @@ -17,22 +17,25 @@ A trivial example of such a map that does not depend on any parameters $a_{n,n+1 ``` In order to solve the discrete Euler-Lagrange equations, the user needs to specify the map $\phi_h$, its gradients with respect to $q_{n}$ and $q_{n+1}$, denoted by $D_1 \phi_h$ and $D_2 \phi_h$, respectively, the gradient with respect to the parameters, denoted by $D_a \phi_h$, and an initial set of parameters $a_{0}$. -The equations of motion, that are solved by this integrator, is computed as: +The equations of motion, that are solved by this integrator, are: ```math \begin{aligned} 0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2}) - + \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n+1}, v_{n+1/2}) \\ - &+ h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} - + h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n} , \\ -0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} , \\ -p_{n+1} &= \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n}, v_{n+1/2}) - + \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n+1}, v_{n+1/2}) , \\ -v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) . + + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + + p_{n} , \\ +0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ +v_{n+1/2} +&= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\ +p_{n+1/2} +&= \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n}, v_{n+1/2}) + + \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n+1}, v_{n+1/2}) , \\ +p_{n+1} +&= h \, D_2 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + + \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n+1}, v_{n+1/2}) . \end{aligned} ``` -The current implementation requires that $D_2 \phi_h $ and $D_a \phi_h$ do not depend on the second and third argument, but only $D_1 \phi_h$ is allowed to depend on all three arguments. -Otherwise it would be necessary to prescribe initial conditions $(q_{-1}, a_{-1,0}, q_{0}, p_{0})$ instead of just $(q_{0}, p_{0})$. - +Given $(q_{n}, p_{n})$, the first four equations are solved for $q_{n+1}$, where $v_{n+1/2}$, $p_{n+1/2}$, and $a_{n,n+1}$ are treated as internal variables similar to the internal stages of a Runge-Kutta method. +The last equation provides an explicit update for $p_{n+1}$. """ struct HPItrapezoidal{ϕT, D₁ϕT, D₂ϕT, DₐϕT, PT} <: HPIMethod ϕ::ϕT From 110b856db500c40246618fccd9c0b4a6793ee6f4 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Wed, 30 Aug 2023 11:58:48 +0900 Subject: [PATCH 09/11] Implement Hamilton-Pontryagin integrators for midpoint Lagrangian. --- docs/src/integrators/hpi.md | 46 ++++++- src/integrators/hpi/hpi_midpoint.jl | 117 +++++++++++++++++- .../hamilton_pontryagin_integrators_tests.jl | 8 +- 3 files changed, 161 insertions(+), 10 deletions(-) diff --git a/docs/src/integrators/hpi.md b/docs/src/integrators/hpi.md index ae98b6278..7715d0a42 100644 --- a/docs/src/integrators/hpi.md +++ b/docs/src/integrators/hpi.md @@ -35,12 +35,11 @@ p_{n+1/2} + \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n+1}, v_{n+1/2}) . \end{aligned} ``` - Upon defining the momentum ```math p_{n} = h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} -+ \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n-1/2}) ++ \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n-1/2}) , ``` we can rewrite the equations of motion as ```math @@ -64,4 +63,47 @@ Given $(q_{n}, p_{n})$, the first four equations can be solved for $q_{n+1}$, wh ## Midpoint Discrete Lagrangian +Similarly to the previous section, we can consider a discrete Lagrangian of the form +```math +L_d (q_{n}, q_{n+1}) = h \, L \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) , +``` +where $q_{n}$ approximates the solution $q(t_n)$ and $v_{n+1/2}$ is the velocity, which is assumed to be constant in the interval $[t_{n}, t_{n+1}]$. +The discrete Hamilton-Pontryagin action reads +```math +A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_{d}^{\alpha,\beta} (q_{n}, q_{n+1}, v_{n+1/2}) + \left< p_{n+1/2} , \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) - v_{n+1/2} \right> \bigg] , +``` +where $\phi_h$ is a map that computes the velocity $v_{n+1/2}$ as a function of $q_{n}$, $q_{n+1}$ and a set of parameters $a_{n,n+1}$. +The equations of motion, that are solved by this integrator, is computed as: +```math +\begin{aligned} +0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n-1} + q_{n}}{2}, v_{n-1/2} \bigg) + + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) \\ + &+ h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + + h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} , \\ +0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ +v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\ +p_{n+1/2} &= \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) . +\end{aligned} +``` +Upon defining the momentum +```math +p_{n} += h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} ++ \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n-1} + q_{n}}{2}, v_{n-1/2} \bigg) , +``` +we can rewrite the equations of motion as +```math +\begin{aligned} +0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) + + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + + p_{n} , \\ +0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} , \\ +v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\ +p_{n+1/2} &= \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) , \\ +p_{n+1} +&= h \, D_2 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) . +\end{aligned} +``` +Given $(q_{n}, p_{n})$, the first four equations can be solved for $q_{n+1}$, where $v_{n+1/2}$, $p_{n+1/2}$, and $a_{n,n+1}$ are treated as internal variables similar to the internal stages of a Runge-Kutta method, and the last equation provides an explicit update for $p_{n+1}$. diff --git a/src/integrators/hpi/hpi_midpoint.jl b/src/integrators/hpi/hpi_midpoint.jl index 8195df2d2..7d0877300 100644 --- a/src/integrators/hpi/hpi_midpoint.jl +++ b/src/integrators/hpi/hpi_midpoint.jl @@ -24,14 +24,33 @@ The equations of motion, that are solved by this integrator, is computed as: + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) \\ &+ h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} , \\ -0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} , \\ +0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ p_{n+1/2} &= \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) , \\ v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) . \end{aligned} ``` -The current implementation requires that $D_2 \phi_h $ and $D_a \phi_h$ do not depend on the second and third argument, but only $D_1 \phi_h$ is allowed to depend on all three arguments. -Otherwise it would be necessary to prescribe initial conditions $(q_{-1}, a_{-1,0}, q_{0}, p_{0})$ instead of just $(q_{0}, p_{0})$. - +Upon defining the momentum +```math +p_{n} += h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} ++ \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n-1} + q_{n}}{2}, v_{n-1/2} \bigg) , +``` +we can rewrite the equations of motion as +```math +\begin{aligned} +0 &= p_{n} + + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) + + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ +0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ +v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\ +p_{n+1/2} &= \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) , \\ +p_{n+1} +&= h \, D_2 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) . +\end{aligned} +``` +Given $(q_{n}, p_{n})$, the first four equations are solved for $q_{n+1}$, where $v_{n+1/2}$, $p_{n+1/2}$, and $a_{n,n+1}$ are treated as internal variables similar to the internal stages of a Runge-Kutta method. +The last equation provides an explicit update for $p_{n+1}$. """ struct HPImidpoint{ϕT, D₁ϕT, D₂ϕT, DₐϕT, PT} <: HPIMethod ϕ::ϕT @@ -50,3 +69,93 @@ issymplectic(method::HPImidpoint) = true const HPImidpointIntegrator{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:HPImidpoint} + +function Cache{ST}(problem::Union{IODEProblem,LODEProblem}, method::HPImidpoint; kwargs...) where {ST} + IntegratorCacheHPI{ST, ndims(problem), nparams(method)}(; kwargs...) +end + +@inline CacheType(ST, problem::Union{IODEProblem,LODEProblem}, method::HPImidpoint) = IntegratorCacheHPI{ST, ndims(problem), nparams(method)} + +function Base.show(io::IO, int::HPImidpointIntegrator) + print(io, "\nHamilton-Pontryagin Integrator using midpoint quadrature with:\n") + print(io, " Timestep: $(timestep(int))\n") +end + + +function components!( + x::Vector{ST}, + solstep::SolutionStepPODE{DT,TT}, + problem::Union{IODEProblem,LODEProblem}, + method::HPImidpoint, + caches::CacheDict) where {ST,DT,TT} + + # get cache and dimension + cache = caches[ST] + D = ndims(problem) + A = nparams(method) + + # set some local variables for convenience and clarity + local t̃ = solstep.t̄ + timestep(problem) / 2 + + # copy x to q + cache.q .= x[1:D] + cache.a .= x[D+1:D+A] + + # compute q̃ + cache.q̃ .= (cache.q .+ cache.q̄) ./ 2 + + # compute v + method.ϕ(cache.ṽ, cache.q̄, cache.q, cache.a, timestep(problem)) + + # compute Θ̃ = ϑ(q̃,ṽ) and f̃ = f(q̃,ṽ) + functions(problem).ϑ(cache.θ̃, t̃, cache.q̃, cache.ṽ) + functions(problem).f(cache.f̃, t̃, cache.q̃, cache.ṽ) + + # compute derivatives of ϕ + method.D₁ϕ(cache.D₁ϕ, cache.q̄, cache.q, cache.a, timestep(problem)) + method.D₂ϕ(cache.D₂ϕ, cache.q̄, cache.q, cache.a, timestep(problem)) + method.Dₐϕ(cache.Dₐϕ, cache.q̄, cache.q, cache.a, timestep(problem)) + + # compute p + cache.p .= timestep(problem) .* cache.f̃ ./ 2 + for i in 1:D + for j in 1:D + cache.p[i] += timestep(problem) * cache.D₂ϕ[i,j] * cache.θ̃[j] + end + end +end + + +function residual!( + b::Vector{ST}, + x::Vector{ST}, + solstep::SolutionStepPODE, + problem::Union{IODEProblem,LODEProblem}, + method::HPImidpoint, + caches::CacheDict) where {ST} + + # get cache for internal stages + cache = caches[ST] + D = ndims(problem) + A = nparams(method) + + # copy previous solution from solstep to cache + reset!(cache, current(solstep)...) + + # compute stages from nonlinear solver solution x + components!(x, solstep, problem, method, caches) + + # compute b + for i in 1:D + b[i] = cache.p̄[i] + timestep(problem) * cache.f̃[i] / 2 + for j in 1:D + b[i] += timestep(problem) * cache.D₁ϕ[i,j] * cache.θ̃[j] + end + end + for i in 1:A + b[D+i] = 0 + for j in 1:D + b[D+i] += cache.Dₐϕ[i,j] * cache.θ̃[j] + end + end +end diff --git a/test/integrators/hamilton_pontryagin_integrators_tests.jl b/test/integrators/hamilton_pontryagin_integrators_tests.jl index 17d0c61d3..2e926765f 100644 --- a/test/integrators/hamilton_pontryagin_integrators_tests.jl +++ b/test/integrators/hamilton_pontryagin_integrators_tests.jl @@ -30,11 +30,11 @@ pref = exact_solution(podeproblem()) end - # sol = integrate(lode, HPImidpoint(ϕ, D₁ϕ, D₂ϕ, Dₐϕ, Float64[])) - # @test relative_maximum_error(sol.q, pref.q) < 4E-4 + sol = integrate(lode, HPImidpoint(ϕ, D₁ϕ, D₂ϕ, Dₐϕ, Float64[])) + @test relative_maximum_error(sol.q, pref.q) < 4E-4 - # ref = integrate(lode, PMVImidpoint()) - # @test relative_maximum_error(sol.q, ref.q) < 8*eps() + ref = integrate(lode, PMVImidpoint()) + @test relative_maximum_error(sol.q, ref.q) < 8*eps() sol = integrate(lode, HPItrapezoidal(ϕ, D₁ϕ, D₂ϕ, Dₐϕ, Float64[])) From e213518795ced2606595dcbb45cbfa2d5ecc87e8 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Wed, 30 Aug 2023 11:59:04 +0900 Subject: [PATCH 10/11] Add missing common files for Hamilton-Pontryagin integrators. --- src/integrators/hpi/hpi_cache.jl | 48 +++++++++++++++++++++++++++++++ src/integrators/hpi/hpi_common.jl | 45 +++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+) create mode 100644 src/integrators/hpi/hpi_cache.jl create mode 100644 src/integrators/hpi/hpi_common.jl diff --git a/src/integrators/hpi/hpi_cache.jl b/src/integrators/hpi/hpi_cache.jl new file mode 100644 index 000000000..4bfd5b6b4 --- /dev/null +++ b/src/integrators/hpi/hpi_cache.jl @@ -0,0 +1,48 @@ +@doc raw""" +Cache for variational integrator in position-momentum form. + +### Fields + +* `q`: internal stages of solution +* `Θ`: implicit function evaluated on solution +* `f`: vector field of implicit function +""" +struct IntegratorCacheHPI{DT,D,A} <: IODEIntegratorCache{DT,D} + x::Vector{DT} + a::Vector{DT} + + q::Vector{DT} + p::Vector{DT} + θ::Vector{DT} + f::Vector{DT} + + q̄::Vector{DT} + p̄::Vector{DT} + θ̄::Vector{DT} + f̄::Vector{DT} + + q̃::Vector{DT} + ṽ::Vector{DT} + θ̃::Vector{DT} + f̃::Vector{DT} + + D₁ϕ::Matrix{DT} + D₂ϕ::Matrix{DT} + Dₐϕ::Matrix{DT} + + function IntegratorCacheHPI{DT,D,A}() where {DT,D,A} + new(zeros(DT, D+A), zeros(DT, A), + zeros(DT,D), zeros(DT,D), zeros(DT,D), zeros(DT,D), + zeros(DT,D), zeros(DT,D), zeros(DT,D), zeros(DT,D), + zeros(DT,D), zeros(DT,D), zeros(DT,D), zeros(DT,D), + zeros(DT,D,D), zeros(DT,D,D), zeros(DT,A,D) + ) + end +end + +function reset!(cache::IntegratorCacheHPI, t, q, p) + copyto!(cache.q̄, q) + copyto!(cache.p̄, p) +end + +nlsolution(cache::IntegratorCacheHPI) = cache.x diff --git a/src/integrators/hpi/hpi_common.jl b/src/integrators/hpi/hpi_common.jl new file mode 100644 index 000000000..d0a20f040 --- /dev/null +++ b/src/integrators/hpi/hpi_common.jl @@ -0,0 +1,45 @@ + +function initial_guess!( + solstep::SolutionStepPODE{DT}, + problem::Union{IODEProblem,LODEProblem}, + method::Union{HPImidpoint,HPItrapezoidal}, + caches::CacheDict, + ::NonlinearSolver, + iguess::Union{InitialGuess,Extrapolation}) where {DT} + + # get cache and dimension + cache = caches[DT] + D = ndims(problem) + A = nparams(method) + + # compute initial guess for solution q(n+1) + initialguess!(solstep.t, cache.q, cache.p, solstep, problem, iguess) + + # copy initial guess to solution vector + cache.x[1:D] .= cache.q + cache.x[D+1:D+A] .= 0 +end + +function integrate_step!( + solstep::SolutionStepPODE{DT,TT}, + problem::Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, + method::Union{HPImidpoint,HPItrapezoidal}, + caches::CacheDict, + solver::NonlinearSolver) where {DT,TT} + + # call nonlinear solver + solve!(caches[DT].x, (b,x) -> residual!(b, x, solstep, problem, method, caches), solver) + + # print solver status + # print_solver_status(int.solver.status, int.solver.params) + + # check if solution contains NaNs or error bounds are violated + # check_solver_status(int.solver.status, int.solver.params) + + # compute vector field at internal stages + components!(caches[DT].x, solstep, problem, method, caches) + + # compute final update + solstep.q .= caches[DT].q + solstep.p .= caches[DT].p +end From e0132c8be6f4522492dff13053ee5804d19e0b18 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Wed, 30 Aug 2023 11:59:16 +0900 Subject: [PATCH 11/11] Minor update in docs of Hamilton-Pontryagin integrators. --- docs/src/integrators/hpi.md | 14 +++++++------- src/integrators/hpi/hpi_trapezoidal.jl | 6 +++--- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/docs/src/integrators/hpi.md b/docs/src/integrators/hpi.md index 7715d0a42..5f8bddbba 100644 --- a/docs/src/integrators/hpi.md +++ b/docs/src/integrators/hpi.md @@ -44,9 +44,9 @@ p_{n} we can rewrite the equations of motion as ```math \begin{aligned} -0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2}) - + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} - + p_{n} , \\ +0 &= p_{n} + + \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2}) + + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ 0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\ @@ -95,10 +95,10 @@ p_{n} we can rewrite the equations of motion as ```math \begin{aligned} -0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) - + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} - + p_{n} , \\ -0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1} , \\ +0 &= p_{n} + + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) + + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ +0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\ p_{n+1/2} &= \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) , \\ p_{n+1} diff --git a/src/integrators/hpi/hpi_trapezoidal.jl b/src/integrators/hpi/hpi_trapezoidal.jl index 1ebe772b0..b39180295 100644 --- a/src/integrators/hpi/hpi_trapezoidal.jl +++ b/src/integrators/hpi/hpi_trapezoidal.jl @@ -20,9 +20,9 @@ In order to solve the discrete Euler-Lagrange equations, the user needs to speci The equations of motion, that are solved by this integrator, are: ```math \begin{aligned} -0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2}) - + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} - + p_{n} , \\ +0 &= p_{n} + + \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2}) + + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ 0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\