Skip to content
Snippets Groups Projects
Commit 17477353 authored by Olivier Cots's avatar Olivier Cots
Browse files

add stabilitiy sys ham

parent c85f6291
No related branches found
No related tags found
No related merge requests found
# Packages needed:
using ForwardDiff: jacobian, gradient, ForwardDiff
using OrdinaryDiffEq: ODEProblem, solve, Tsit5, OrdinaryDiffEq
# --------------------------------------------------------------------------------------------
# Default options for flows
# --------------------------------------------------------------------------------------------
__abstol() = 1e-10
__reltol() = 1e-10
__saveat() = []
__method() = OrdinaryDiffEq.Tsit5()
# --------------------------------------------------------------------------------------------------
# A desription is a tuple of symbols
const DescVarArg = Vararg{Symbol} # or Symbol...
const Description = Tuple{DescVarArg}
# --------------------------------------------------------------------------------------------------
# the description may be given as a tuple or a list of symbols (Vararg{Symbol})
makeDescription(desc::DescVarArg) = Tuple(desc) # create a description from Vararg{Symbol}
makeDescription(desc::Description) = desc
# default is autonomous
isnonautonomous(desc::Description) = :nonautonomous desc
# --------------------------------------------------------------------------------------------------
# Aliases for types
#
const Time = Number
const State = Vector{<:Number} # Vector{de sous-type de Number}
const Adjoint = Vector{<:Number}
const CoTangent = Vector{<:Number}
const Control = Vector{<:Number}
const DState = Vector{<:Number}
const DAdjoint = Vector{<:Number}
const DCoTangent = Vector{<:Number}
# --------------------------------------------------------------------------------------------
# Hamiltonian
# --------------------------------------------------------------------------------------------
struct Hamiltonian f::Function end
function (h::Hamiltonian)(x::State, p::Adjoint, λ...) # https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects
return h.f(x, p, λ...)
end
function (h::Hamiltonian)(t::Time, x::State, p::Adjoint, λ...)
return h.f(t, x, p, λ...)
end
# Flow from a Hamiltonian
# On peut mettre les options d'intégration dans Flow
# ou à l'appel de f par la suite
function Flow(h::Hamiltonian, description...; kwargs_Flow...)
h_(t, x, p, λ...) = isnonautonomous(makeDescription(description...)) ?
h(t, x, p, λ...) : h(x, p, λ...)
function rhs!(dz::DCoTangent, z::CoTangent, λ, t::Time)
n = size(z, 1)÷2
foo = z -> h_(t, z[1:n], z[n+1:2*n], λ...)
dh = ForwardDiff.gradient(foo, z)
dz[1:n] = dh[n+1:2n]
dz[n+1:2n] = -dh[1:n]
end
function f(tspan::Tuple{Time, Time}, x0::State, p0::Adjoint, λ...;
method=__method(), abstol=__abstol(), reltol=__reltol(), saveat=__saveat(),
kwargs...)
z0 = [ x0 ; p0 ]
ode = OrdinaryDiffEq.ODEProblem(rhs!, z0, tspan, λ)
sol = OrdinaryDiffEq.solve(ode, method, abstol=abstol, reltol=reltol, saveat=saveat;
kwargs..., kwargs_Flow...)
return sol
end
function f(t0::Time, x0::State, p0::Adjoint, tf::Time, λ...;
method=__method(), abstol=__abstol(), reltol=__reltol(), saveat=__saveat(),
kwargs...)
sol = f((t0, tf), x0, p0, λ..., method=method, abstol=abstol, reltol=reltol, saveat=saveat;
kwargs..., kwargs_Flow...)
n = size(x0, 1)
return sol[1:n, end], sol[n+1:2*n, end]
end
return f
end;
\ No newline at end of file
td/stability/pendule-inv.png

19.8 KiB

Source diff could not be displayed: it is too large. Options to address this: view the blob.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment