Examples
Unstable manifold of Henon map
Consider the Henon map:
\[\begin{aligned} x'&=1-\alpha x^2+y,\\ y'&=\beta x, \end{aligned}\]
where $\alpha=1.4,\beta=0.3$. This map has a saddle located at $(0.6313544770895048, 0.18940634312685142)$, and its unstable eigenvector is $(-0.9880577559947047, 0.15408397327012555)$.
To compute the unstable manifolds of the saddle numerically, InvariantManifolds.jl just needs a segment of unstable manifold, whose start point is the saddle. It's resonable to choose a short unstable eigenvector as the segment. You don't have to shorten the eigenvector started at the saddle yourself. We provide a function segment
to do this automatically. The segment
can generate equal distributed n
points at one point, with given length and direction.
using InvariantManifolds, StaticArrays
function henonmap(x, p)
y1 = 1 - p[1] * x[1]^2 + x[2]
y2 = p[2] * x[1]
SA[y1, y2]
end
henonmap2(x, p)=henonmap(henonmap(x, p), p)
# For technical reason, we have to use the double iteration of the map, since the eigenvalue is less than -1.
seg = segment(SA[0.6313544770895048, 0.18940634312685142], SA[-0.9880577559947047, 0.15408397327012555], 150, 0.01)
result = generate_curves(henonmap2, SA[1.4, 0.3], seg, 0.001, 8)
This package does not provide any function to plot the manifolds. However, it's simple to plot it using the stardard julia ploting library.
To do this, just run
using GLMakie
function manifold_plot(v::Vector{IterationCurve})
figure = Figure()
axes = Axis(figure[1,1])
for k in eachindex(v)
data = v[k].pcurve.u
lines!(axes, first.(data), last.(data))
end
figure
end
manifold_plot(result)
Unstable manifold of the periodic perturbed system:
Consider the Duffing equation with periodic perturbation:
\[\begin{aligned} \dot{x}&=y,\\ \dot{y}&=x-x^3+\gamma \cos(t). \end{aligned}\]
When $\gamma=0$, the system has a saddle located at $(0,0)$. After the small periodic perturbation, this saddle becomes a saddle periodic orbit, i.e., a saddle of the map $T:x\mapsto \phi(X,2\pi,0)$, where $\phi(X,t,t_0)$ is the solution of system with initial condition $X(t_0)=X\in\mathbb{R}^2$. Fortunately, we can use the solution of the variational equation to get the jacobian matrix of $T$. The map $T$'s saddle's location and unstable direction can also be obtained numerically.
using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, GLMakie
# Duffing Equation
f(x, p, t) = SA[x[2], x[1] - x[1]^3 + p[1]*cos(t)]
# First find the fixed point and its unstable direction by using the Newton iteration.
function timemap(x,p)
prob = ODEProblem{false}(f, x, (0.0, 2pi), p)
solve(prob, Vern9())[end]
end
# Solving the variational equation to get the jacobian matrix.
function jac(x, p)
prob = ODEProblem{false}(f, x, (0.0, 2pi), p)
sol = solve(prob, Vern9())
function df(x, p, t)
SA[0 1; 1-3*(sol(t)[1])^2 0] * x
end
ii = SA[1.0 0.0; 0.0 1.0]
nprob = ODEProblem{false}(df, ii, (0.0, 2pi), p)
solve(nprob, Vern9())[end]
end
# Newton's iteration to get the saddle's location.
function newton(x, p; n=100, atol=1e-8)
xn = x - inv(jac(x, p) - I) * (timemap(x, p) - x)
data = typeof(x)[x, xn]
i = 1
while norm(data[2] - data[1]) > atol && i <= n
data[1] = data[2]
data[2] = data[1] - inv(jac(data[1], p) - I) * (timemap(data[1], p) - data[1])
i = i + 1
end
if norm(data[2] - data[1]) < atol
println("Fixed point found successfully:")
data[2]
else
println("Failed to find a fixed point after $n times iterations. The last point is:")
data[2]
end
end
function manifold_plot(v)
figure = Figure()
axes = Axis(figure[1,1])
for k in eachindex(v)
data = v[k].pcurve.u
lines!(axes, first.(data), last.(data))
end
figure
end
para = [0.1]
fixedpoint = newton(SA[-0.05, 0.0], para)
unstable_direction = eigen(jac(fixedpoint, para)).vectors[:, 2]
seg = segment(fixedpoint, unstable_direction, 150, 0.01)
result = generate_curves(timemap, para, seg, 0.002, 3)
manifold_plot(result)
Lorenz manifold
In this example, we will consider the well known Lorenz's equation. With classical parameters, the fixed point origin has two linear independent stable direction. Hence, there exists a two-dimensional stable manifolds of the origin, the so called Lorenz manifold. The function generate_surface
mainly needs two parameter to compute such stable manifolds. The first parameter is the time-$T$-map of the system, where $T<0$. Moreover, the original vector field has to be rescaled to ensure the uniform extension of the flow. The second parameter is the two linear independent stable direction. See the blew codes for more details.
using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, GLMakie
# First define the rescaled lorenz vector field and its time-T-map
function lorenz(x, p, t)
σ, ρ, β = p
v = SA[σ*(x[2]-x[1]),
ρ*x[1]-x[2]-x[1]*x[3],
x[1]*x[2]-β*x[3]
]
v / sqrt(1 + norm(v)^2)
end
function lorenz_map(x, p)
prob = ODEProblem{false}(lorenz, x, (0.0, -2.0), p)
sol = solve(prob, Tsit5())
sol[end]
end
function eigenv(p)
σ, ρ, β = p
(SA[0.0, 0.0, 1.0], SA[-(-1 + σ + sqrt(1 - 2 * σ + 4 * ρ * σ + σ^2))/(2*ρ), 1, 0])
end
second(x) = x[2]
function manifold_plot(annulus)
fig = Figure()
axes = LScene(fig[1, 1], show_axis=false,scenekw = (backgroundcolor=:white, clear=true))
for i in eachindex(annulus)
points = annulus[i].outer.pcurve.u
scatter!(axes, first.(points), second.(points), last.(points), fxaa=true)
end
fig
end
para = [10, 28, 9 / 3]
lorenz_manifold = generate_surface(lorenz_map, para, SA[0.0, 0.0, 0.0], eigenv(para)..., 120, 1, 1)
manifold_plot(lorenz_manifold)
Unstable manifold of the piecewise smooth ODE's time-$T$-map
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 & \text{if } x < -d,\\ k_2 & \text{if } -d<x<d,\\ -k_3 & \text{if } x > d. \end{cases}\]
where $k_1,k_2,k_3,d>0$ are positive parameters. We are going to compute the time-1-map' unstable manifolds of saddle located near $(0,0)$ when $\epsilon$ is small.
First, we define a piecewise smooth vector field:
using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, GLMakie
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), 0)
The parameters pass to PiecewiseV
are vector fields, their definition domains, and the hyper surfaces separating these domains. See PiecewiseV
for more details.
Then we set the information needs when computing the time-1-map:
setup = setmap(vectorfield, (0.0, 1.0), Tsit5(), 2, Float64)
See setmap
for the meaning of parameters of this function. For small $\epsilon$, the amplitude of the saddle periodic orbit is small so that the orbit is in domain dom1
. So we can still using the method of variational equation to compute the saddle's location.
function timemap(x,p)
prob = ODEProblem{false}(f1, x, (0.0, 1.0), p)
solve(prob, Vern9())[end]
end
function jac(x,p)
prob = ODEProblem{false}(f1, x, (0.0, 1.0), p)
sol = solve(prob, Vern9())
function df(x, p, t)
SA[0 1; p[1] 0] * x
end
ii = SA[1.0 0.0; 0.0 1.0]
nprob = ODEProblem{false}(df, ii, (0.0, 1.0), p)
solve(nprob, Vern9())[end]
end
function newton(x,p; n=100, atol=1e-8)
xn = x - inv(jac(x,p) - I) * (timemap(x,p) - x)
data = [x, xn]
i = 1
while norm(data[2] - data[1]) > atol && i <= n
data[1] = data[2]
data[2] = data[1] - inv(jac(data[1],p) - I) * (timemap(data[1],p) - data[1])
end
if norm(data[2] - data[1]) < atol
println("Fixed point found successfully:")
data[2]
else
println("Failed to find a fixed point after $n times iterations. The last point is:")
data[2]
end
end
Then we define the parameters and try to find the saddle:
para = [2, 5, 5, 0.6, 2]
fixedpoint= newton(SA[0.0, 0.0],para)
unstable_direction = eigen(jac(fixedpoint,para)).vectors[:,2]
As before, we compute the manifold and plot it using GLMakie
. Noting that the data structures of result
is slightly different from examples before.
seg = segment(fixedpoint, unstable_direction, 150, 0.05)
result = generate_curves(setup, para, seg, 0.01, 9)
function manifold_plot(result)
fig = Figure()
axes = Axis(fig[1,1])
for k in eachindex(result)
for j in eachindex(result[k])
data=result[k][j].pcurve.u
lines!(axes,first.(data),last.(data))
end
end
fig
end
manifold_plot(result)
Unstable manifold of the impact inverted pendulum's time-$T$-map
Consider an inverted pendulum with two side walls and a small periodic perturbation:
\[\begin{aligned} \dot{x}&= y,\\ \dot{y}&= sin(x) - \epsilon \cos(2\pi t), \end{aligned}\]
When $x=\xi$ or $x=-\xi$, we have $y\rightarrow - ry$.
The procedure for computing such system's manifold is similar with before, we just give codes here:
using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, GLMakie
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(), 2, Float64)
# First find the fixed point and its unstable direction
function timemap(x,p)
prob = ODEProblem{false}(f, x, (0.0, 1.0), p)
sol = solve(prob, Vern9())[end]
end
function jac(x, p)
prob = ODEProblem{false}(f, x, (0.0, 1.0), p)
sol = solve(prob, Vern9())
function df(x, p, t)
SA[0 1; cos(sol(t)[1]) 0] * x
end
ii = SA[1.0 0.0; 0.0 1.0]
nprob = ODEProblem{false}(df, ii, (0.0, 1.0), p)
solve(nprob, Vern9())[end]
end
function newton(x, p; n=100, atol=1e-12)
xn = x - inv(jac(x, p) - I) * (timemap(x, p) - x)
data = typeof(x)[x, xn]
i = 1
while norm(data[2] - data[1]) > atol && i <= n
data[1] = data[2]
data[2] = data[1] - inv(jac(data[1], p) - I) * (timemap(data[1], p) - data[1])
i = i + 1
end
if norm(data[2] - data[1]) < atol
println("Fixed point found successfully:")
data[2]
else
println("Failed to find a fixed point after $n times iterations. The last point is:")
data[2]
end
end
para = [0.2, pi / 4, 0.98, 2]
fixedpoint = newton(SA[0.0, 0.0], para)
unstable_direction = eigen(jac(fixedpoint, para)).vectors[:, 2]
seg = segment(fixedpoint, unstable_direction, 150, 0.01)
result = generate_curves(setup, para, seg, 0.001, 11)
function manifold_plot(result)
fig = Figure()
axes = Axis(fig[1,1])
for k in eachindex(result)
for j in eachindex(result[k])
data=result[k][j].pcurve.u
lines!(axes,first.(data),last.(data))
end
end
fig
end
manifold_plot(result)
Unstable manifold of the Filippov system's time-$T$-map
I don't found any existing simple examples, but I think the codes should just work.