One-Dimensional Non-Smooth Manifolds
Perhaps the most notable feature of this package is its ability to compute non-smooth manifolds. Currently, it supports the computation of non-smooth manifolds for two types of systems:
- One-dimensional manifold computation for time-periodic non-smooth differential equations
- Two-dimensional manifold computation for non-smooth autonomous systems
The manifolds in both cases are invariant manifolds of saddle points. The former requires taking time-periodic mappings, while the latter requires taking fixed-step time-T-mappings. The non-smooth factors in these two types of systems can be diverse, including piecewise smooth functions, collisions, and combinations thereof. Users don't need to solve these three types of non-smooth systems themselves, as we provide three encapsulated structures:
We use the Callback functionality from OrdinaryDiffEq.jl to compute time mappings. The following three examples will demonstrate the methods for computing invariant manifolds of these three types of non-smooth systems.
The computation of non-smooth manifolds heavily depends on the ODE solving algorithm and precision. When the solution fails or performs poorly, try changing the algorithm, increasing the ODE solving precision, or reducing the values of parameters ϵ
, d
, amax
in NSOneDManifoldProblem
.
Piecewise Smooth Systems
Consider a simple piecewise smooth system:
\[\begin{aligned} \dot{x}&=y,\\ \dot{y}&=f(x) + \epsilon \sin(2\pi t), \end{aligned}\]
where
\[f(x) = \begin{cases} -k_1 x& \text{if } x < -d,\\ k_2 x & \text{if } -d<x<d,\\ -k_3 x& \text{if } x > d. \end{cases}\]
\[k_1,k_2,k_3,d>0\]
are all positive constants. We will compute the invariant manifold of the time-periodic mapping. Note that when the periodic perturbation is small, the saddle point should be close to the origin.
First, let's load the required packages
julia> using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, CairoMakie
Next, define the piecewise smooth vector field:
f1(x, p, t) = SA[x[2], p[1]*x[1]+p[4]*sin(2pi * t)]
f2(x, p, t) = SA[x[2], -p[2]*x[1]+p[4]*sin(2pi * t)]
f3(x, p, t) = SA[x[2], -p[3]*x[1]+p[4]*sin(2pi * t)]
hyper1(x, p, t) = x[1] - p[5]
hyper2(x, p, t) = x[1] + p[5]
dom1(x, p, t) = -p[5] < x[1] < p[5]
dom2(x, p, t) = x[1] > p[5]
dom3(x, p, t) = x[1] < -p[5]
vectorfield = PiecewiseV((f1, f2, f3), (dom1, dom2, dom3), (hyper1, hyper2))
PiecewiseV:
fs: (Main.f1, Main.f2, Main.f3)
regions: (Main.dom1, Main.dom2, Main.dom3)
hypers: (Main.hyper1, Main.hyper2)
n: 0
The parameters passed to the PiecewiseV
structure are: vector fields, their domains, and the hyperplanes that divide these domains. For more details, refer to PiecewiseV
.
Next, we'll encapsulate the key information for solving the time-periodic mapping in another structure NSSetUp
:
julia> setup = setmap(vectorfield, (0.0, 1.0), Tsit5(), abstol=1e-8)
NSSetUp: f: PiecewiseV{Tuple{typeof(Main.f1), typeof(Main.f2), typeof(Main.f3)}, Tuple{typeof(Main.dom1), typeof(Main.dom2), typeof(Main.dom3)}, Tuple{typeof(Main.hyper1), typeof(Main.hyper2)}}((Main.f1, Main.f2, Main.f3), (Main.dom1, Main.dom2, Main.dom3), (Main.hyper1, Main.hyper2), 0) timespan: (0.0, 1.0)
The function setmap
is used to encapsulate the time mapping computation information. Now we have defined everything needed to solve the time-periodic mapping.
Next, to generate the local manifold, we also need to locate the saddle point and its unstable eigenvector. Let's set the parameters:
julia> para = [2, 5, 5, 0.6, 2]
5-element Vector{Float64}: 2.0 5.0 5.0 0.6 2.0
Since the perturbation is small, the saddle-type periodic orbit should still be in dom1
. Therefore, we can use findsaddle
to calculate the position of the saddle point:
function df1(x, p, t)
SA[0 1; p[1] 0]
end
initialguess = SA[0.0, 0.0]
saddle = findsaddle(f1, df1, (0.0,1.0), initialguess, para, abstol=1e-10)
Saddle{2,Float64,Float64}:
saddle: [1.3120846754625683e-11, -0.0908885006640104]
unstable_directions: StaticArraysCore.SVector{2, Float64}[[0.5773502691896258, 0.816496580927726]]
unstable_eigen_values: [4.113250379341709]
Next, create the problem, generate the local manifold, and perform the extension
julia> prob = NSOneDManifoldProblem(setup, para, ϵ = 1e-3)
NSOneDManifoldProblem: f: NSSetUp generated by non-smooth vector field: PiecewiseV{Tuple{typeof(Main.f1), typeof(Main.f2), typeof(Main.f3)}, Tuple{typeof(Main.dom1), typeof(Main.dom2), typeof(Main.dom3)}, Tuple{typeof(Main.hyper1), typeof(Main.hyper2)}}((Main.f1, Main.f2, Main.f3), (Main.dom1, Main.dom2, Main.dom3), (Main.hyper1, Main.hyper2), 0) para: [2.0, 5.0, 5.0, 0.6, 2.0] amax: 0.5 d: 0.001 ϵ: 0.001 dsmin: 1.0e-6
julia> segment = gen_segment(saddle)
50-element Vector{StaticArraysCore.SVector{2, Float64}}: [1.3120846754625683e-11, -0.0908885006640104] [0.00011782659866974998, -0.09072186870871903] [0.00023565318421865322, -0.09055523675342765] [0.0003534797697675564, -0.09038860479813628] [0.0004713063553164597, -0.0902219728428449] [0.000589132940865363, -0.09005534088755354] [0.0007069595264142661, -0.08988870893226217] [0.0008247861119631695, -0.0897220769769708] [0.0009426126975120726, -0.08955544502167942] [0.001060439283060976, -0.08938881306638805] ⋮ [0.0048308900206258795, -0.08405659049706411] [0.004948716606174783, -0.08388995854177275] [0.005066543191723685, -0.08372332658648138] [0.005184369777272589, -0.08355669463119] [0.005302196362821493, -0.08339006267589863] [0.005420022948370395, -0.08322343072060726] [0.005537849533919299, -0.08305679876531588] [0.005655676119468201, -0.08289016681002451] [0.005773502705017105, -0.08272353485473313]
julia> manifold = growmanifold(prob, segment, 8)
Non-smooth one-dimensional manifold Curves number: 30 Points number: 158521 Total arc length: 91.63468743370944 Flaw points number: 11 Distance failed points number: 0 Curvature failed points number: 11 Max dα in Flaw Points: 1.246076029203152e-6
Note that the data type of manifold.data
is Vector{Vector{S}}
, where S
is an interpolation function. So we need to use the following function to plot the results:
using CairoMakie
function manifold_plot(data)
fig = Figure()
axes = Axis(fig[1,1])
for k in eachindex(data)
for j in eachindex(data[k])
points=Point2f.(data[k][j].u)
lines!(axes, points)
end
end
fig
end
manifold_plot(manifold.data)

Full codes without comments:
using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, CairoMakie
f1(x, p, t) = SA[x[2], p[1]*x[1]+p[4]*sin(2pi * t)]
f2(x, p, t) = SA[x[2], -p[2]*x[1]+p[4]*sin(2pi * t)]
f3(x, p, t) = SA[x[2], -p[3]*x[1]+p[4]*sin(2pi * t)]
hyper1(x, p, t) = x[1] - p[5]
hyper2(x, p, t) = x[1] + p[5]
dom1(x, p, t) = -p[5] < x[1] < p[5]
dom2(x, p, t) = x[1] > p[5]
dom3(x, p, t) = x[1] < -p[5]
vectorfield = PiecewiseV((f1, f2, f3), (dom1, dom2, dom3), (hyper1, hyper2))
setup = setmap(vectorfield, (0.0, 1.0), Tsit5(), abstol=1e-8)
para = [2, 5, 5, 0.6, 2]
function df1(x, p, t)
SA[0 1; p[1] 0]
end
initialguess = SA[0.0, 0.0]
saddle = findsaddle(f1, df1, (0.0,1.0), initialguess, para, abstol=1e-10)
prob = NSOneDManifoldProblem(setup, para, ϵ = 1e-3)
segment = gen_segment(saddle)
manifold = growmanifold(prob, segment, 8)
using CairoMakie
function manifold_plot(data)
fig = Figure()
axes = Axis(fig[1,1])
for k in eachindex(data)
for j in eachindex(data[k])
points=Point2f.(data[k][j].u)
lines!(axes,points)
end
end
fig
end
manifold_plot(manifold.data)
Impact Systems
Consider the following forced inverted pendulum equation:
\[\begin{aligned} \dot{x}&= y,\\ \dot{y}&= \sin(x) - \epsilon \cos(2\pi t), \end{aligned}\]
Assume there are walls on both sides of the inverted pendulum: when $x=\xi$ or $x=-\xi$, we have $y\rightarrow - ry$. Similarly, we need to first construct the non-smooth vector field
f(x, p, t) = SA[x[2], sin(x[1])-p[1]*cos(2 * pi * t)]
hyper1(x, p, t) = x[1] + p[2]
hyper2(x, p, t) = x[1] - p[2]
rule1(x, p, t) = SA[x[1], -p[3]*x[2]]
rule2(x, p, t) = SA[x[1], -p[3]*x[2]]
vectorfield = BilliardV(f, (hyper1, hyper2), (rule1, rule2))
BilliardV{typeof(Main.f), Tuple{typeof(Main.hyper1), typeof(Main.hyper2)}, Tuple{typeof(Main.rule1), typeof(Main.rule2)}}(Main.f, (Main.hyper1, Main.hyper2), (Main.rule1, Main.rule2))
Next, encapsulate the information for solving the time-periodic mapping:
setup = setmap(vectorfield, (0.0, 1.0), Vern9(), abstol=1e-10)
NSSetUp:
f: BilliardV{typeof(Main.f), Tuple{typeof(Main.hyper1), typeof(Main.hyper2)}, Tuple{typeof(Main.rule1), typeof(Main.rule2)}}(Main.f, (Main.hyper1, Main.hyper2), (Main.rule1, Main.rule2))
timespan: (0.0, 1.0)
Find the saddle point:
function df(x, p, t)
SA[0 1; cos(x[1]) 0]
end
para = [0.2, pi / 4, 0.98]
initialguess = SA[0.0, 0.0]
saddle = findsaddle(f, df, (0.0,1.0), initialguess, para, abstol=1e-10)
Saddle{2,Float64,Float64}:
saddle: [0.004940905039667357, 3.526850741022572e-11]
unstable_directions: StaticArraysCore.SVector{2, Float64}[[0.707107886757561, 0.7071056756138054]]
unstable_eigen_values: [2.7182735331462027]
Next, create the problem, generate the local manifold, and perform the extension
julia> prob = NSOneDManifoldProblem(setup, para)
NSOneDManifoldProblem: f: NSSetUp generated by non-smooth vector field: BilliardV{typeof(Main.f), Tuple{typeof(Main.hyper1), typeof(Main.hyper2)}, Tuple{typeof(Main.rule1), typeof(Main.rule2)}}(Main.f, (Main.hyper1, Main.hyper2), (Main.rule1, Main.rule2)) para: [0.2, 0.7853981633974483, 0.98] amax: 0.5 d: 0.001 ϵ: 1.0e-5 dsmin: 1.0e-6
julia> segment = gen_segment(saddle)
50-element Vector{StaticArraysCore.SVector{2, Float64}}: [0.004940905039667357, 3.526850741022572e-11] [0.005085212771658696, 0.00014430731600601873] [0.0052295205036500345, 0.00028861459674353005] [0.005373828235641373, 0.00043292187748104134] [0.005518135967632713, 0.0005772291582185527] [0.005662443699624052, 0.000721536438956064] [0.005806751431615391, 0.0008658437196935753] [0.00595105916360673, 0.001010151000431087] [0.006095366895598068, 0.001154458281168598] [0.006239674627589407, 0.0012987655619061092] ⋮ [0.010857522051312256, 0.005916598545506473] [0.011001829783303593, 0.006060905826243983] [0.011146137515294932, 0.006205213106981494] [0.011290445247286271, 0.006349520387719005] [0.011434752979277612, 0.006493827668456517] [0.011579060711268949, 0.0066381349491940285] [0.01172336844326029, 0.00678244222993154] [0.011867676175251626, 0.006926749510669051] [0.012011983907242967, 0.007071056791406562]
julia> manifold = growmanifold(prob, segment, 11)
Non-smooth one-dimensional manifold Curves number: 38 Points number: 96016 Total arc length: 63.81546496564522 Flaw points number: 3 Distance failed points number: 0 Curvature failed points number: 3 Max dα in Flaw Points: 2.4211133345433726e-7
Finally, plot the results:
manifold_plot(manifold.data)

Full codes without comments:
using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, CairoMakie
f(x, p, t) = SA[x[2], sin(x[1])-p[1]*cos(2 * pi * t)]
hyper1(x, p, t) = x[1] + p[2]
hyper2(x, p, t) = x[1] - p[2]
rule1(x, p, t) = SA[x[1], -p[3]*x[2]]
rule2(x, p, t) = SA[x[1], -p[3]*x[2]]
vectorfield = BilliardV(f, (hyper1, hyper2), (rule1, rule2))
setup = setmap(vectorfield, (0.0, 1.0), Vern9(), abstol=1e-10)
function df(x, p, t)
SA[0 1; cos(x[1]) 0]
end
para = [0.2, pi / 4, 0.98]
initialguess = SA[0.0, 0.0]
saddle = findsaddle(f, df, (0.0,1.0), initialguess, para, abstol=1e-10)
prob = NSOneDManifoldProblem(setup, para)
segment = gen_segment(saddle)
manifold = growmanifold(prob, segment, 11)
function manifold_plot(data)
fig = Figure()
axes = Axis(fig[1,1])
for k in eachindex(data)
for j in eachindex(data[k])
points=Point2f.(data[k][j].u)
lines!(axes,points)
end
end
fig
end
manifold_plot(manifold.data)
Combined Piecewise Smooth and Impact ODE Systems
Now consider an ODE system with both piecewise smooth functions and impacts:
\[\begin{aligned} \dot{x}&=y,\\ \dot{y}&=f(x) + \epsilon \sin(2\pi t), \end{aligned}\]
where
\[f(x) = \begin{cases} -k_1 x& \text{if } x < -d,\\ k_2 x & \text{if } -d<x<d \end{cases}\]
\[k_1,k_2,d>0\]
are all positive constants. When $x=d$, we have $\dot{y}->-r\dot{y}$. We will compute the invariant manifold of the time-periodic mapping. Note that when the periodic perturbation is small, the saddle point should be close to the origin. First, let's load the required packages
julia> using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, CairoMakie
Next, define the non-smooth vector field:
f1(x, p, t) = SA[x[2], p[1]*x[1]+p[3]*sin(2pi * t)]
f2(x, p, t) = SA[x[2], -p[2]*x[1]+p[3]*sin(2pi * t)]
hyper1(x, p, t) = x[1] - p[4]
hyper2(x, p, t) = x[1] + p[4]
dom1(x, p, t) = -p[4] < x[1]
dom2(x, p, t) = x[1] < -p[4]
impact_rule(x, p, t) = SA[x[1], -p[5]*x[2]]
id(x,p,t) = x
vectorfield = PiecewiseImpactV((f1, f2), (dom1, dom2), (hyper1, hyper2), (impact_rule, id), [1])
PiecewiseImpactV{Tuple{typeof(Main.f1), typeof(Main.f2)}, Tuple{typeof(Main.hyper1), typeof(Main.hyper2)}, Tuple{typeof(Main.impact_rule), typeof(Main.id)}, Tuple{typeof(Main.dom1), typeof(Main.dom2)}}((Main.f1, Main.f2), (Main.dom1, Main.dom2), (Main.hyper1, Main.hyper2), (Main.impact_rule, Main.id), [1], 0)
The parameters passed to the PiecewiseImpactV
structure are: vector fields, their domains, the hyperplanes that divide these domains, the rules that act on the hyperplanes, and a list of rules that have impact effects. For more details, refer to PiecewiseImpactV
.
Next, we'll encapsulate the key information for solving the time-periodic mapping in another structure NSSetUp
:
julia> setup = setmap(vectorfield, (0.0, 1.0), Tsit5(), abstol=1e-8, reltol=1e-8)
NSSetUp: f: PiecewiseImpactV{Tuple{typeof(Main.f1), typeof(Main.f2)}, Tuple{typeof(Main.hyper1), typeof(Main.hyper2)}, Tuple{typeof(Main.impact_rule), typeof(Main.id)}, Tuple{typeof(Main.dom1), typeof(Main.dom2)}}((Main.f1, Main.f2), (Main.dom1, Main.dom2), (Main.hyper1, Main.hyper2), (Main.impact_rule, Main.id), [1], 0) timespan: (0.0, 1.0)
The function setmap
is used to encapsulate the time mapping computation information. Now we have defined everything needed to solve the time-periodic mapping.
Next, to generate the local manifold, we also need to locate the saddle point and its unstable eigenvector. Let's set the parameters:
julia> para = [2, 5, 0.5, 2, 0.98]
5-element Vector{Float64}: 2.0 5.0 0.5 2.0 0.98
Since the perturbation is small, the saddle-type periodic orbit should still be in dom1
. Therefore, we can use findsaddle
to calculate the position of the saddle point:
function df1(x, p, t)
SA[0 1; p[1] 0]
end
initialguess = SA[0.0, 0.0]
saddle = findsaddle(f1, df1, (0.0, 1.0), initialguess, para, abstol=1e-10)
Saddle{2,Float64,Float64}:
saddle: [7.410052858050141e-12, -0.0757404172152893]
unstable_directions: StaticArraysCore.SVector{2, Float64}[[0.5773502691896258, 0.816496580927726]]
unstable_eigen_values: [4.113250379341709]
Next, create the problem, generate the local manifold, and perform the extension
julia> prob = NSOneDManifoldProblem(setup, para)
NSOneDManifoldProblem: f: NSSetUp generated by non-smooth vector field: PiecewiseImpactV{Tuple{typeof(Main.f1), typeof(Main.f2)}, Tuple{typeof(Main.hyper1), typeof(Main.hyper2)}, Tuple{typeof(Main.impact_rule), typeof(Main.id)}, Tuple{typeof(Main.dom1), typeof(Main.dom2)}}((Main.f1, Main.f2), (Main.dom1, Main.dom2), (Main.hyper1, Main.hyper2), (Main.impact_rule, Main.id), [1], 0) para: [2.0, 5.0, 0.5, 2.0, 0.98] amax: 0.5 d: 0.001 ϵ: 1.0e-5 dsmin: 1.0e-6
julia> segment = gen_segment(saddle)
50-element Vector{StaticArraysCore.SVector{2, Float64}}: [7.410052858050141e-12, -0.0757404172152893] [0.0001178265929589561, -0.07557378525999793] [0.00023565317850785932, -0.07540715330470656] [0.00035347976405676254, -0.07524052134941518] [0.0004713063496056658, -0.07507388939412381] [0.000589132935154569, -0.07490725743883245] [0.0007069595207034722, -0.07474062548354107] [0.0008247861062523756, -0.0745739935282497] [0.0009426126918012787, -0.07440736157295832] [0.001060439277350182, -0.07424072961766695] ⋮ [0.004830890014915086, -0.06890850704834302] [0.004948716600463989, -0.06874187509305166] [0.0050665431860128915, -0.06857524313776028] [0.005184369771561795, -0.06840861118246891] [0.005302196357110699, -0.06824197922717753] [0.0054200229426596015, -0.06807534727188616] [0.005537849528208505, -0.06790871531659479] [0.005655676113757407, -0.06774208336130341] [0.0057735026993063114, -0.06757545140601204]
julia> manifold = growmanifold(prob, segment, 9)
Non-smooth one-dimensional manifold Curves number: 52 Points number: 292201 Total arc length: 175.48287334895852 Flaw points number: 2 Distance failed points number: 0 Curvature failed points number: 2 Max dα in Flaw Points: 1.4116867281496803e-7
Note that the data type of manifold.data
is Vector{Vector{S}}
, where S
is an interpolation function. So we need to use the following function to plot the results:
using CairoMakie
function manifold_plot(data)
fig = Figure()
axes = Axis(fig[1,1])
for k in eachindex(data)
for j in eachindex(data[k])
points=Point2f.(data[k][j].u)
lines!(axes,points)
end
end
fig
end
manifold_plot(manifold.data)

Full codes without comments:
using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, CairoMakie
f1(x, p, t) = SA[x[2], p[1]*x[1]+p[3]*sin(2pi * t)]
f2(x, p, t) = SA[x[2], -p[2]*x[1]+p[3]*sin(2pi * t)]
hyper1(x, p, t) = x[1] - p[4]
hyper2(x, p, t) = x[1] + p[4]
dom1(x, p, t) = -p[4] < x[1]
dom2(x, p, t) = x[1] < -p[4]
impact_rule(x, p, t) = SA[x[1], -p[5]*x[2]]
id(x,p,t) = x
vectorfield = PiecewiseImpactV((f1, f2), (dom1, dom2), (hyper1, hyper2), (impact_rule, id), [1])
setup = setmap(vectorfield, (0.0, 1.0), Tsit5(), abstol=1e-8, reltol=1e-8)
para = [2, 5, 0.5, 2, 0.98]
initialguess = SA[0.0, 0.0]
function df1(x, p, t)
SA[0 1; p[1] 0]
end
saddle = findsaddle(f1, df1, (0.0,1.0), initialguess, para, abstol=1e-10)
segment = gen_segment(saddle)
prob = NSOneDManifoldProblem(setup, para)
manifold = growmanifold(prob, segment, 9)
function manifold_plot(data)
fig = Figure()
axes = Axis(fig[1,1])
for k in eachindex(data)
for j in eachindex(data[k])
points=Point2f.(data[k][j].u)
lines!(axes,points)
end
end
fig
end
manifold_plot(manifold.data)