Two-Dimensional Smooth Manifolds

Essentially, we haven't introduced new algorithms; the core functions for computing two-dimensional manifolds are the same as those for one-dimensional manifolds. We simply represent the two-dimensional manifold as circles of one-dimensional manifolds that are close enough to each other.

Autonomous Vector Field: Lorenz Manifold

First, let's load the required packages and define the Lorenz vector field:

julia> using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, CairoMakie
julia> 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(0.1 + norm(v)^2)
       endlorenz (generic function with 1 method)

It's worth noting that we performed an approximate normalization of the vector field to keep its magnitude within a small range. This ensures uniform expansion of the manifold. Under the classical parameters:

julia> para = [10.0, 28.0, 8/3]3-element Vector{Float64}:
 10.0
 28.0
  2.6666666666666665

the Jacobian matrix at the origin equilibrium point has two stable directions, which are:

function eigenv(p)
    σ, ρ, β = p
    [SA[0.0, 0.0, 1.0], SA[-(-1 + σ + sqrt(1 - 2 * σ + 4 * ρ * σ + σ^2))/(2*ρ), 1, 0]]
end
eigenv(para)

2-element Vector{StaticArraysCore.SVector{3, Float64}}:
 [0.0, 0.0, 1.0]
 [-0.7795615518272664, 1.0, 0.0]

Then we can create a Saddle structure to store this saddle point:

julia> saddle = Saddle(SA[0, 0, 0.0], eigenv(para), [1.0, 1.0])Saddle{3,Float64,Float64}:
saddle: [0.0, 0.0, 0.0]
unstable_directions: StaticArraysCore.SVector{3, Float64}[[0.0, 0.0, 1.0], [-0.7795615518272664, 1.0, 0.0]]
unstable_eigen_values: [1.0, 1.0]

The magnitude of the eigenvalues can be specified arbitrarily, as it won't affect the computation results. Since we're computing the stable manifold, we need to evolve the flow backward. Let's define the following mapping:

julia> function lorenz_map(x, p)
           prob = ODEProblem{false}(lorenz, x, (0.0, -1.0), p)
           sol = solve(prob, Vern9(), abstol = 1e-10)
           sol[end]
       endlorenz_map (generic function with 1 method)

Now we can create the problem:

julia> prob = VTwoDManifoldProblem(lorenz_map, para, d=1.0, amax=1.0, dsmin=1e-3)VTwoDManifoldProblem:
f:Main.lorenz_map
para: [10.0, 28.0, 2.6666666666666665]
amax: 1.0
d: 1.0
dsmin: 0.001

For the meaning of these keyword arguments, please refer to VTwoDManifoldProblem.

Similar to one-dimensional manifolds, we need a local manifold to start the extension. The corresponding function for creating a local manifold is gen_disk:

julia> disk = gen_disk(saddle, r=1.0)2-element Vector{Vector{StaticArraysCore.SVector{3, Float64}}}:
 [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
 [[0.0, 0.0, 1.0], [-0.00010127432293440568, 0.00012991190073063764, 0.9999999864332046], [-0.00020254864312087533, 0.00025982379793629897, 0.9999999457328191], [-0.00030382295781147296, 0.0003897356880920077, 0.9999998778988444], [-0.00040509726425826294, 0.0005196475676727878, 0.9999997829312824], [-0.0005063715597133095, 0.0006495594331536636, 0.9999996608301356], [-0.0006076458414286776, 0.0007794712810096597, 0.9999995115954075], [-0.000708920106656432, 0.000909383107715801, 0.9999993352271019], [-0.0008101943526486386, 0.0010392949097471134, 0.9999991317252238], [-0.0009114685766573632, 0.001169206683578623, 0.9999989010897786]  …  [0.0009114685766575245, -0.0011692066835788302, 0.9999989010897786], [0.0008101943526485436, -0.0010392949097469916, 0.9999991317252238], [0.0007089201066560809, -0.0009093831077153506, 0.9999993352271019], [0.000607645841428499, -0.0007794712810094307, 0.9999995115954075], [0.0005063715597133038, -0.0006495594331536562, 0.9999996608301356], [0.0004050972642584299, -0.000519647567673002, 0.9999997829312824], [0.0003038229578118127, -0.0003897356880924434, 0.9999998778988444], [0.00020254864312095882, -0.00025982379793640603, 0.9999999457328191], [0.00010127432293423296, -0.0001299119007304161, 0.9999999864332046], [0.0, 0.0, 1.0]]

For detailed information about the gen_disk function, please refer to gen_disk. Now we can proceed with the extension:

julia> manifold = growmanifold(prob, disk, 200)Two-dimensional manifold
Circles number: 202
Points number: 8028260
Flaw points number: 0

We can also define a plotting function to visualize the results:

using CairoMakie
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 = Point3f.(annulus[i].u)
        lines!(axes, points, fxaa=true)
    end
    fig
end
manifold_plot(manifold.data)

Full codes without comments:

using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, CairoMakie
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(0.1 + norm(v)^2)
end
function eigenv(p)
    σ, ρ, β = p
    [SA[0.0, 0.0, 1.0], SA[-(-1 + σ + sqrt(1 - 2 * σ + 4 * ρ * σ + σ^2))/(2*ρ), 1, 0]]
end
para = [10.0, 28.0, 8/3]
saddle = Saddle(SA[0.0, 0.0, 0.0], eigenv(para), [1.0, 1.0])
prob = VTwoDManifoldProblem(lorenz, para, d=1.0, amax=1.0, dsmin=1e-3)
disk = gen_disk(saddle, r=1.0)
manifold = growmanifold(prob, disk, 200)
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 = Point3f.(annulus[i].u)
        lines!(axes, points, fxaa=true)
    end
    fig
end
manifold_plot(manifold.data)

Nonlinear Mapping

Consider the following nonlinear mapping:

\[f(X)=\varphi\circ\Lambda\circ\varphi^{-1}(X)\]

where $\varphi(x,y,z)=(x,y,z-\alpha x^2-\beta y^2)$ is a nonlinear mapping, and $\Lambda$ is a diagonal matrix whose diagonal elements can be used to control the Jacobian matrix of mapping $f$ near the origin. Here's the code for computing its invariant manifold:

using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, CairoMakie
const Λ = SDiagonal(SA[2.1, 6.3, 0.6])
φ(x, p)= SA[x[1],x[2],x[3]-p[1]*x[1]^2-p[2]*x[2]^2]
iφ(x, p)= SA[x[1],x[2],x[3]+p[1]*x[1]^2+p[2]*x[2]^2]
f(x,p) = φ(Λ*iφ(x, p),p)

para = [1.2,-1.2]
saddle = Saddle(SA[0.0, 0.0, 0.0], [SA[1.0, 0.0, 0.0], SA[0.0, 1.0, 0.0]], [2.1, 6.3])
prob = TwoDManifoldProblem(f, para, dcircle=0.05, d = 0.02, dsmin=1e-3)

disk = gen_disk(saddle, times=4, r= 0.05)
manifold = growmanifold(prob, disk, 3)
function manifold_plot(data)
    fig = Figure()
    axes = Axis3(fig[1,1])
    for k in eachindex(data)
        for j in eachindex(data[k])
            points=data[k][j].u
            scatter!(axes,Point3f.(points),fxaa=true)
        end
    end
    fig
end
manifold_plot(manifold.data)