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.

Warning

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)