Getting Started: One-Dimensional Smooth Manifolds
Nonlinear Mapping
Consider the following Henon map:
\[\begin{aligned} x'&=1-\alpha x^2+y,\\ y'&=\beta x, \end{aligned}\]
where $\alpha,\beta$ are parameters. This mapping has fixed points:
\[\begin{aligned} (x_1,y_1)&=(\frac{-\sqrt{4 \alpha +\beta ^2-2 \beta +1}+\beta -1}{2 \alpha },\frac{1}{2} \left(\frac{\beta ^2}{\alpha }-\frac{\beta \sqrt{4 \alpha +\beta ^2-2 \beta +1}}{\alpha }-\frac{\beta }{\alpha }\right)),\\ (x_2,y_2)&=(\frac{\sqrt{4 \alpha +\beta ^2-2 \beta +1}+\beta -1}{2 \alpha },\frac{1}{2} \left(\frac{\beta ^2}{\alpha }+\frac{\beta \sqrt{4 \alpha +\beta ^2-2 \beta +1}}{\alpha }-\frac{\beta }{\alpha }\right)), \end{aligned}\]
Let's calculate the eigenvalues of these two fixed points under the classical parameters $\alpha=1.4,\beta=0.3$:
using StaticArrays, LinearAlgebra
function fixedpoints(p)
a , b = p
x1 = (-sqrt(4 * a + b^2 - 2 * b + 1) + b - 1) / (2 * a)
y1 = (1 / 2) * (b^2 / a - b * sqrt(4 * a + b^2 - 2 * b + 1) / a - b / a)
x2 = (sqrt(4 * a + b^2 - 2 * b + 1) + b - 1) / (2 * a)
y2 = (1 / 2) * (b^2 / a + b * sqrt(4 * a + b^2 - 2 * b + 1) / a - b / a)
return SA[x1, y1], SA[x2, y2]
end
function jacobian(x, p)
a, b = p
J = @SMatrix [-2 * a * x[1] 1.0; b 0.0]
return J
endjacobian (generic function with 1 method)julia> eigen(jacobian(fixedpoints([1.4, 0.3])[1], [1.4, 0.3]))LinearAlgebra.Eigen{Float64, Float64, StaticArraysCore.SMatrix{2, 2, Float64, 4}, StaticArraysCore.SVector{2, Float64}} values: 2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2): -0.0920295620408392 3.259822097891452 vectors: 2×2 StaticArraysCore.SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2): -0.293276 0.995792 0.956028 0.0916423
julia> eigen(jacobian(fixedpoints([1.4, 0.3])[2], [1.4, 0.3]))LinearAlgebra.Eigen{Float64, Float64, StaticArraysCore.SMatrix{2, 2, Float64, 4}, StaticArraysCore.SVector{2, Float64}} values: 2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2): -1.9237388581534067 0.15594632230279395 vectors: 2×2 StaticArraysCore.SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2): -0.988058 -0.461228 0.154084 -0.887282
As we can see, under the classical parameters, both fixed points are unstable. Next, we'll use the InvariantManifolds.jl package to compute one branch of the unstable manifold of the second fixed point.
The InvariantManifolds.jl package has an interface similar to many Julia packages. First, we need to load the package in Julia and define the Henon map:
julia> using InvariantManifoldsjulia> function henonmap(x, p) y1 = 1 - p[1] * x[1]^2 + x[2] y2 = p[2] * x[1] SA[y1, y2] endhenonmap (generic function with 1 method)
Since the unstable eigenvalue at the saddle point is:
julia> eigen(jacobian(fixedpoints([1.4, 0.3])[2], [1.4, 0.3])).values[1]-1.9237388581534067
We need to iterate this mapping twice to ensure the manifold doesn't reverse during extension:
julia> henonmap2(x, p)=henonmap(henonmap(x, p), p)henonmap2 (generic function with 1 method)
Now let's define a problem for computing a one-dimensional manifold of the smooth mapping:
julia> para = [1.4, 0.3]2-element Vector{Float64}: 1.4 0.3julia> prob = OneDManifoldProblem(henonmap2, para)OneDManifoldProblem: f:Main.henonmap2 para: [1.4, 0.3] amax: 0.5 d: 0.001 dsmin: 1.0e-6
To compute the manifold, we need a small local manifold segment starting at the saddle point. Usually, an unstable eigenvector starting at the saddle point with a very small length will suffice. InvariantManifolds.jl provides a function gen_segment to generate such a local manifold:
saddle = fixedpoints(para)[2]
unstable_direction = eigen(jacobian(fixedpoints([1.4, 0.3])[2], [1.4, 0.3])).vectors[:,1]
segment = gen_segment(saddle, unstable_direction)50-element Vector{StaticArraysCore.SVector{2, Float64}}:
[0.6313544770895047, 0.1894063431268514]
[0.6311528326495057, 0.18943778883568202]
[0.6309511882095068, 0.18946923454451267]
[0.6307495437695079, 0.1895006802533433]
[0.6305478993295089, 0.18953212596217395]
[0.63034625488951, 0.18956357167100457]
[0.6301446104495111, 0.18959501737983522]
[0.6299429660095122, 0.18962646308866585]
[0.6297413215695133, 0.1896579087974965]
[0.6295396771295144, 0.18968935450632712]
⋮
[0.623087055049549, 0.19069561718890754]
[0.62288541060955, 0.1907270628977382]
[0.6226837661695511, 0.19075850860656882]
[0.6224821217295522, 0.19078995431539947]
[0.6222804772895533, 0.1908214000242301]
[0.6220788328495543, 0.19085284573306074]
[0.6218771884095554, 0.19088429144189137]
[0.6216755439695565, 0.19091573715072202]
[0.6214738995295576, 0.19094718285955264]Under default keyword arguments, this function will generate a local manifold starting at the saddle point, with a length of 150 units and a step size of 0.01. Now let's use this local manifold to compute the smooth manifold:
julia> manifold = growmanifold(prob, segment, 8)One-dimensional manifold Curves number: 10 Points number: 264346 Total arc length: 174.42126112443742 Flaw points number: 39 Distance failed points number: 0 Curvature failed points number: 39 Max dα in Flaw Points: 1.7498998144196122e-5
This package doesn't provide plotting functionality. However, since the computation results are stored in manifold.data, and manifold.data is actually a vector whose elements are interpolation functions from the DataInterpolations.jl package:
julia> manifold.data10-element Vector{DataInterpolations.QuadraticInterpolation{Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Float64}, Nothing, Vector{StaticArraysCore.SVector{2, Float64}}, StaticArraysCore.SVector{2, Float64}, (2,)}}: DataInterpolations.QuadraticInterpolation{Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Float64}, Nothing, Vector{StaticArraysCore.SVector{2, Float64}}, StaticArraysCore.SVector{2, Float64}, (2,)}(StaticArraysCore.SVector{2, Float64}[[0.6313544770895047, 0.1894063431268514], [0.6311528326495057, 0.18943778883568202], [0.6309511882095068, 0.18946923454451267], [0.6307495437695079, 0.1895006802533433], [0.6305478993295089, 0.18953212596217395], [0.63034625488951, 0.18956357167100457], [0.6301446104495111, 0.18959501737983522], [0.6299429660095122, 0.18962646308866585], [0.6297413215695133, 0.1896579087974965], [0.6295396771295144, 0.18968935450632712] … [0.6232886994895479, 0.19066417148007692], [0.623087055049549, 0.19069561718890754], [0.62288541060955, 0.1907270628977382], [0.6226837661695511, 0.19075850860656882], [0.6224821217295522, 0.19078995431539947], [0.6222804772895533, 0.1908214000242301], [0.6220788328495543, 0.19085284573306074], [0.6218771884095554, 0.19088429144189137], [0.6216755439695565, 0.19091573715072202], [0.6214738995295576, 0.19094718285955264]], [0.0, 0.0002040816326530667, 0.0004081632653061377, 0.0006122448979592045, 0.0008163265306122755, 0.0010204081632653422, 0.0012244897959184133, 0.00142857142857148, 0.0016326530612244413, 0.001836734693877508 … 0.008163265306122424, 0.00836734693877549, 0.008571428571428561, 0.008775510204081629, 0.0089795918367347, 0.009183673469387766, 0.009387755102040837, 0.009591836734693904, 0.009795918367346975, 0.010000000000000042], nothing, DataInterpolations.QuadraticParameterCache{Vector{StaticArraysCore.SVector{2, Float64}}}(StaticArraysCore.SVector{2, Float64}[], StaticArraysCore.SVector{2, Float64}[]), :Forward, DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{Vector{Float64}}([0.0, 0.0002040816326530667, 0.0004081632653061377, 0.0006122448979592045, 0.0008163265306122755, 0.0010204081632653422, 0.0012244897959184133, 0.00142857142857148, 0.0016326530612244413, 0.001836734693877508 … 0.008163265306122424, 0.00836734693877549, 0.008571428571428561, 0.008775510204081629, 0.0089795918367347, 0.009183673469387766, 0.009387755102040837, 0.009591836734693904, 0.009795918367346975, 0.010000000000000042], Base.RefValue{Int64}(1), true), false, false) DataInterpolations.QuadraticInterpolation{Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Float64}, Nothing, Vector{StaticArraysCore.SVector{2, Float64}}, StaticArraysCore.SVector{2, Float64}, (2,)}(StaticArraysCore.SVector{2, Float64}[[0.6214738995295576, 0.19094718285955264], [0.6201363288165692, 0.19114810131287846], [0.6193867223329736, 0.19126394528685345], [0.6186369017026461, 0.19137975510602506], [0.6178868672947623, 0.19149553077039336], [0.6171366194783905, 0.1916112722799583], [0.6163861586224891, 0.19172697963471985], [0.6156354850959078, 0.19184265283467805], [0.6148845992673879, 0.19195829187983285], [0.614133501505561, 0.19207389677018435] … [0.6013327704354634, 0.19403395422124683], [0.6005779207860675, 0.194148944325138], [0.5998228661866863, 0.1942639002742258], [0.5990676070037746, 0.19437882206851023], [0.5983121436036772, 0.19449370970799132], [0.5975564763526315, 0.19460856319266906], [0.5968006056167654, 0.19472338252254343], [0.5960445317620974, 0.19483816769761447], [0.5952882551545376, 0.19495291871788215], [0.5945317761598876, 0.1950676355833465]], [0.0, 0.0013525766658977183, 0.002111081581211276, 0.00286979291754507, 0.003628710315862905, 0.004387833417215384, 0.0051471618627421615, 0.0059066952936710415, 0.006666433351317443, 0.007426375677086039 … 0.02037630118750546, 0.021139859120596308, 0.021903614525341418, 0.022667567044922364, 0.023431716322611085, 0.02419606200176783, 0.024960603725842856, 0.025725341138376503, 0.026490273882998225, 0.027255401603427105], nothing, DataInterpolations.QuadraticParameterCache{Vector{StaticArraysCore.SVector{2, Float64}}}(StaticArraysCore.SVector{2, Float64}[], StaticArraysCore.SVector{2, Float64}[]), :Forward, DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{Vector{Float64}}([0.0, 0.0013525766658977183, 0.002111081581211276, 0.00286979291754507, 0.003628710315862905, 0.004387833417215384, 0.0051471618627421615, 0.0059066952936710415, 0.006666433351317443, 0.007426375677086039 … 0.02037630118750546, 0.021139859120596308, 0.021903614525341418, 0.022667567044922364, 0.023431716322611085, 0.02419606200176783, 0.024960603725842856, 0.025725341138376503, 0.026490273882998225, 0.027255401603427105], Base.RefValue{Int64}(35), false), false, true) DataInterpolations.QuadraticInterpolation{Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Float64}, Nothing, Vector{StaticArraysCore.SVector{2, Float64}}, StaticArraysCore.SVector{2, Float64}, (2,)}(StaticArraysCore.SVector{2, Float64}[[0.5945317761598876, 0.1950676355833465], [0.5939066367276442, 0.1951623537504024], [0.5932812144166023, 0.1952570733893607], [0.5926555091695513, 0.19535179450749376], [0.5920295209291975, 0.1954465171120734], [0.5914032496381673, 0.19554124121037045], [0.5907766952390041, 0.1956359668096553], [0.5901498576741722, 0.1957306939171972], [0.5895227368860531, 0.1958254225402649], [0.5888187542042903, 0.195931717664183] … [0.4985086193049846, 0.20913517525928219], [0.49777989497957087, 0.2092384350950152], [0.4970508538251686, 0.20934169039181894], [0.49632168223887285, 0.20944491477359053], [0.49559238052039334, 0.2095481082400438], [0.49486294896935534, 0.2096512707908924], [0.4941332019089832, 0.20975442870941446], [0.4934033255492317, 0.20985755569599662], [0.49267332018969856, 0.2099606517503523], [0.49194318612989485, 0.21006371687219508]], [0.0, 0.0006322743399158579, 0.001264828588303605, 0.001897662805462797, 0.00253077705177201, 0.0031641713876851426, 0.003797845873735015, 0.0044318005705294315, 0.005066035538754482, 0.005777997805725589 … 0.09704852001192218, 0.0977845239087361, 0.09852084086571197, 0.09925728263208132, 0.09999384891429618, 0.10073053941888226, 0.10146754167686767, 0.10220466763346231, 0.102941916995192, 0.10367928946865838], nothing, DataInterpolations.QuadraticParameterCache{Vector{StaticArraysCore.SVector{2, Float64}}}(StaticArraysCore.SVector{2, Float64}[], StaticArraysCore.SVector{2, Float64}[]), :Forward, DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{Vector{Float64}}([0.0, 0.0006322743399158579, 0.001264828588303605, 0.001897662805462797, 0.00253077705177201, 0.0031641713876851426, 0.003797845873735015, 0.0044318005705294315, 0.005066035538754482, 0.005777997805725589 … 0.09704852001192218, 0.0977845239087361, 0.09852084086571197, 0.09925728263208132, 0.09999384891429618, 0.10073053941888226, 0.10146754167686767, 0.10220466763346231, 0.102941916995192, 0.10367928946865838], Base.RefValue{Int64}(144), true), false, false) DataInterpolations.QuadraticInterpolation{Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Float64}, Nothing, Vector{StaticArraysCore.SVector{2, Float64}}, StaticArraysCore.SVector{2, Float64}, (2,)}(StaticArraysCore.SVector{2, Float64}[[0.49194318612989485, 0.21006371687219508], [0.49133974412388703, 0.21014886103070418], [0.4907362143635332, 0.21023398410535546], [0.49013259701732814, 0.21031908609599798], [0.4895288922537252, 0.21040416700248107], [0.48892483292856626, 0.21048926447508395], [0.48832068644173054, 0.21057434084519858], [0.48771645296176935, 0.21065939611267423], [0.487112132657195, 0.21074443027736006], [0.48650745817259866, 0.21082948096028176] … [0.0914486599132987, 0.2606262198742126], [0.09071778890201432, 0.26070959321032167], [0.08998681198305208, 0.2607929518005496], [0.08925585172308409, 0.2608762816866813], [0.08852490834694027, 0.26095958286852133], [0.08779398207937703, 0.26104285534587435], [0.08706295109357329, 0.2611261130168524], [0.0863319376715399, 0.2612093419733651], [0.08560094203791957, 0.26129254221521714], [0.0848699644172827, 0.2613757137422129]], [0.0, 0.0006094192172412676, 0.0012189223833288891, 0.001828509333198196, 0.002438179901821013, 0.003048203875247368, 0.003658311215235842, 0.004268501756650052, 0.00487877533438766, 0.005489401939595204 … 0.40368539822883975, 0.4044210092322307, 0.40515672379007245, 0.4058924185444205, 0.4066280932728328, 0.4073637477529377, 0.4080995046027066, 0.4088352407526919, 0.4095709559806274, 0.41030665006431627], nothing, DataInterpolations.QuadraticParameterCache{Vector{StaticArraysCore.SVector{2, Float64}}}(StaticArraysCore.SVector{2, Float64}[], StaticArraysCore.SVector{2, Float64}[]), :Forward, DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{Vector{Float64}}([0.0, 0.0006094192172412676, 0.0012189223833288891, 0.001828509333198196, 0.002438179901821013, 0.003048203875247368, 0.003658311215235842, 0.004268501756650052, 0.00487877533438766, 0.005489401939595204 … 0.40368539822883975, 0.4044210092322307, 0.40515672379007245, 0.4058924185444205, 0.4066280932728328, 0.4073637477529377, 0.4080995046027066, 0.4088352407526919, 0.4095709559806274, 0.41030665006431627], Base.RefValue{Int64}(576), true), false, false) DataInterpolations.QuadraticInterpolation{Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Float64}, Nothing, Vector{StaticArraysCore.SVector{2, Float64}}, StaticArraysCore.SVector{2, Float64}, (2,)}(StaticArraysCore.SVector{2, Float64}[[0.0848699644172827, 0.2613757137422129], [0.0842658442218916, 0.261444431167246], [0.08366173660123466, 0.26151312897926254], [0.083057641681961, 0.2615818071781525], [0.0824535595906854, 0.261650465763806], [0.08184940724898274, 0.2617191141891766], [0.08124526798805473, 0.2617877429961905], [0.08064114193445951, 0.26185635218473785], [0.0800370292147185, 0.26192494175470893], [0.07943284692126251, 0.2619935211296817] … [-1.1614732345562178, 0.3750360898261848], [-1.1620437498454024, 0.37507558762223364], [-1.16261306123063, 0.3751149767725769], [-1.1631810888334355, 0.37515425177474443], [-1.1637478162911987, 0.3751934115270798], [-1.1643132592059073, 0.375232457137918], [-1.1648774011194916, 0.37527138750028305], [-1.1654402577351426, 0.37531020372789786], [-1.1660018125104026, 0.3753489047094255], [-1.1665620812197832, 0.37538749156138373]], [0.0, 0.0006080158673772994, 0.001216017024211429, 0.0018240033451883183, 0.0024319747050267973, 0.003040014718751775, 0.003648039520669254, 0.004256048985550897, 0.004864042988203632, 0.005472104970568422 … 1.251598951108218, 1.2521708320146014, 1.252741504388661, 1.2533108881637585, 1.2538789669413815, 1.2544457563639362, 1.2550112399373283, 1.2555754334053058, 1.256138320189294, 1.2566999161043921], nothing, DataInterpolations.QuadraticParameterCache{Vector{StaticArraysCore.SVector{2, Float64}}}(StaticArraysCore.SVector{2, Float64}[], StaticArraysCore.SVector{2, Float64}[]), :Forward, DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{Vector{Float64}}([0.0, 0.0006080158673772994, 0.001216017024211429, 0.0018240033451883183, 0.0024319747050267973, 0.003040014718751775, 0.003648039520669254, 0.004256048985550897, 0.004864042988203632, 0.005472104970568422 … 1.251598951108218, 1.2521708320146014, 1.252741504388661, 1.2533108881637585, 1.2538789669413815, 1.2544457563639362, 1.2550112399373283, 1.2555754334053058, 1.256138320189294, 1.2566999161043921], Base.RefValue{Int64}(1922), false), false, false) DataInterpolations.QuadraticInterpolation{Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Float64}, Nothing, Vector{StaticArraysCore.SVector{2, Float64}}, StaticArraysCore.SVector{2, Float64}, (2,)}(StaticArraysCore.SVector{2, Float64}[[-1.1665620812197832, 0.37538749156138373], [-1.1674853308769213, 0.37545102169915395], [-1.1684050417388405, 0.3755142385019023], [-1.169321212038843, 0.3755771419799027], [-1.1702338400225136, 0.3756397321436884], [-1.1711430489226031, 0.37570201756060884], [-1.1720487110633844, 0.37576398959926777], [-1.1729508247265181, 0.37582564827097403], [-1.1738493882059564, 0.3758869935872952], [-1.174744522585078, 0.3759480339274975] … [0.26337041013990464, -0.15773752025608406], [0.2626662190197014, -0.1578723565085279], [0.26196301081164003, -0.15800689762930156], [0.26125915508906156, -0.15814145589040238], [0.2605546516365122, -0.15827603129186932], [0.25984950023851383, -0.15841062383374077], [0.2591453200321457, -0.1585449247959122], [0.2584404944202113, -0.158679242820039], [0.2577350231886236, -0.15881357790616013], [0.25702890612326884, -0.1589479300543143]], [0.0, 0.0009254328759076409, 0.0018473137975466123, 0.0027656410083767985, 0.003680412764161521, 0.004591752600426186, 0.005499532555492402, 0.006403750921269816, 0.0073044060019831734, 0.008201619176104907 … 2.658716959571874, 2.6594339434969005, 2.6601499065525562, 2.6608665088159054, 2.6615837505138944, 2.6623016318734702, 2.6630185046037935, 2.6637360145165085, 2.664454161837022, 2.6651729467907455], nothing, DataInterpolations.QuadraticParameterCache{Vector{StaticArraysCore.SVector{2, Float64}}}(StaticArraysCore.SVector{2, Float64}[], StaticArraysCore.SVector{2, Float64}[]), :Forward, DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{Vector{Float64}}([0.0, 0.0009254328759076409, 0.0018473137975466123, 0.0027656410083767985, 0.003680412764161521, 0.004591752600426186, 0.005499532555492402, 0.006403750921269816, 0.0073044060019831734, 0.008201619176104907 … 2.658716959571874, 2.6594339434969005, 2.6601499065525562, 2.6608665088159054, 2.6615837505138944, 2.6623016318734702, 2.6630185046037935, 2.6637360145165085, 2.664454161837022, 2.6651729467907455], Base.RefValue{Int64}(4031), false), false, false) DataInterpolations.QuadraticInterpolation{Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Float64}, Nothing, Vector{StaticArraysCore.SVector{2, Float64}}, StaticArraysCore.SVector{2, Float64}, (2,)}(StaticArraysCore.SVector{2, Float64}[[0.25702890612326884, -0.1589479300543143], [0.25644662877791435, -0.15905863980790858], [0.2558639127048849, -0.15916936115033628], [0.25528075778431397, -0.15928009408161914], [0.2546971638963211, -0.1593908386017791], [0.2541131309210128, -0.15950159471083813], [0.2535286587384841, -0.1596123624088179], [0.25294374722881524, -0.15972314169574056], [0.2523583962720775, -0.1598339325716274], [0.2517748550100404, -0.1599443097240967] … [0.29712097571348495, 0.22406084680956365], [0.2966190326498853, 0.2241176724847115], [0.29611835873013315, 0.22417432287919337], [0.2956182490799543, 0.2242308778655098], [0.2951182603972722, 0.224287387590453], [0.2946188382041869, 0.2243438017083997], [0.29411954004940954, 0.22440017021391448], [0.29362081060708445, 0.224456442913375], [0.2931222083540081, 0.2245126696385835], [0.29262417689196474, 0.22456880037887758]], [0.0, 0.000592708660687505, 0.0011858504965378029, 0.0017794256336995432, 0.002373434198324772, 0.0029678763165678773, 0.003562752114584198, 0.004158061718533459, 0.004753805254573546, 0.005347693729194098 … 7.169624754040204, 7.170129903519519, 7.170633772193992, 7.171137069450769, 7.171640241421196, 7.172142839749536, 7.172645309704635, 7.173147203790421, 7.173648966335257, 7.174150150929186], nothing, DataInterpolations.QuadraticParameterCache{Vector{StaticArraysCore.SVector{2, Float64}}}(StaticArraysCore.SVector{2, Float64}[], StaticArraysCore.SVector{2, Float64}[]), :Forward, DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{Vector{Float64}}([0.0, 0.000592708660687505, 0.0011858504965378029, 0.0017794256336995432, 0.002373434198324772, 0.0029678763165678773, 0.003562752114584198, 0.004158061718533459, 0.004753805254573546, 0.005347693729194098 … 7.169624754040204, 7.170129903519519, 7.170633772193992, 7.171137069450769, 7.171640241421196, 7.172142839749536, 7.172645309704635, 7.173147203790421, 7.173648966335257, 7.174150150929186], Base.RefValue{Int64}(11033), true), false, true) DataInterpolations.QuadraticInterpolation{Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Float64}, Nothing, Vector{StaticArraysCore.SVector{2, Float64}}, StaticArraysCore.SVector{2, Float64}, (2,)}(StaticArraysCore.SVector{2, Float64}[[0.29262417689196474, 0.22456880037887758], [0.2918040740999408, 0.2246611612247733], [0.2909849306822663, 0.22475332808854487], [0.2901667525985359, 0.22484530035164216], [0.28934954581691175, 0.22493707739479787], [0.2885333163140996, 0.22502865859802706], [0.2877180700753288, 0.22512004334062602], [0.2869038130943214, 0.22521123100117282], [0.2860905513732776, 0.2253022209575255], [0.28528140667316315, 0.2253926644843337] … [-0.6135234092751721, 0.33077674568223525], [-0.6143216979587726, 0.33084703138611027], [-0.6151193757659539, 0.3309172446584382], [-0.6159166462219929, 0.33098740345184635], [-0.6167125959506898, 0.33105742745386246], [-0.617508138501074, 0.3311273970969567], [-0.6183030651887255, 0.3311972940655287], [-0.619097584152727, 0.33126713669883623], [-0.6198907805782264, 0.3313368446425278], [-0.6206835695013369, 0.3314064983749493]], [0.0, 0.0008252872926078707, 0.0016495995329009478, 0.0024729307383559085, 0.0032952749178713583, 0.004116626071791433, 0.004936978191926354, 0.005756325261583215, 0.006574661255580894, 0.007388845004660185 … 14.873532275413377, 14.874333652292588, 14.875134414296982, 14.875934765742514, 14.876733789729307, 14.877532403340778, 14.878330397090727, 14.879127979925036, 14.87992423350703, 14.8807200764055], nothing, DataInterpolations.QuadraticParameterCache{Vector{StaticArraysCore.SVector{2, Float64}}}(StaticArraysCore.SVector{2, Float64}[], StaticArraysCore.SVector{2, Float64}[]), :Forward, DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{Vector{Float64}}([0.0, 0.0008252872926078707, 0.0016495995329009478, 0.0024729307383559085, 0.0032952749178713583, 0.004116626071791433, 0.004936978191926354, 0.005756325261583215, 0.006574661255580894, 0.007388845004660185 … 14.873532275413377, 14.874333652292588, 14.875134414296982, 14.875934765742514, 14.876733789729307, 14.877532403340778, 14.878330397090727, 14.879127979925036, 14.87992423350703, 14.8807200764055], Base.RefValue{Int64}(23597), false), false, false) DataInterpolations.QuadraticInterpolation{Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Float64}, Nothing, Vector{StaticArraysCore.SVector{2, Float64}}, StaticArraysCore.SVector{2, Float64}, (2,)}(StaticArraysCore.SVector{2, Float64}[[-0.6206835695013369, 0.3314064983749493], [-0.6213359955517115, 0.3314638061539772], [-0.6219881448479169, 0.33152107717631357], [-0.6226400171902914, 0.33157831144164956], [-0.6232916123792942, 0.33163550894967603], [-0.6239421609346247, 0.3316926021946916], [-0.6245924325482655, 0.331749658764159], [-0.6252424270217224, 0.3318066786577696], [-0.6258921441566281, 0.3318636618752149], [-0.6265408110740576, 0.33192054067118065] … [-0.07050334164246772, 0.23849981878217855], [-0.0698356118048438, 0.23840177715407887], [-0.06916904511194633, 0.23830386546853458], [-0.06850252898907142, 0.2382059203572876], [-0.06783606374141546, 0.23810794182028866], [-0.06716964967427438, 0.23800992985748753], [-0.06650363076940127, 0.23791193504629723], [-0.0658376636023805, 0.23781390684367568], [-0.06517174847835105, 0.23771584524957295], [-0.06450588570254889, 0.23761775026393803]], [0.0, 0.0006549381136752951, 0.0013095973185602737, 0.0019639774159665047, 0.002618078207325188, 0.0032711272571438453, 0.003923897219401843, 0.004576387896570421, 0.005228599091246255, 0.005879754961371057 … 42.903706532422284, 42.90438142152198, 42.90505514094002, 42.905728815183906, 42.90640244395411, 42.907076026951025, 42.907749216496896, 42.90842235971795, 42.90909545631477, 42.909768505987834], nothing, DataInterpolations.QuadraticParameterCache{Vector{StaticArraysCore.SVector{2, Float64}}}(StaticArraysCore.SVector{2, Float64}[], StaticArraysCore.SVector{2, Float64}[]), :Forward, DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{Vector{Float64}}([0.0, 0.0006549381136752951, 0.0013095973185602737, 0.0019639774159665047, 0.002618078207325188, 0.0032711272571438453, 0.003923897219401843, 0.004576387896570421, 0.005228599091246255, 0.005879754961371057 … 42.903706532422284, 42.90438142152198, 42.90505514094002, 42.905728815183906, 42.90640244395411, 42.907076026951025, 42.907749216496896, 42.90842235971795, 42.90909545631477, 42.909768505987834], Base.RefValue{Int64}(62550), true), false, true) DataInterpolations.QuadraticInterpolation{Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Float64}, Nothing, Vector{StaticArraysCore.SVector{2, Float64}}, StaticArraysCore.SVector{2, Float64}, (2,)}(StaticArraysCore.SVector{2, Float64}[[-0.06450588570254889, 0.23761775026393803], [-0.0639579540826353, 0.23753699804723477], [-0.06341005829279417, 0.23745622321628498], [-0.06286219850330682, 0.23737542577106058], [-0.062314374884505, 0.23729460571153385], [-0.06176682087667673, 0.23721379746970853], [-0.061219303349430904, 0.23713296663281733], [-0.06067182247302361, 0.2370521132008323], [-0.060124378417760055, 0.23697123717372573], [-0.059577204727296745, 0.23689037304655933] … [-1.144520872853663, 0.3694622427378967], [-1.1444348680282457, 0.3694721878222378], [-1.144344581330131, 0.369481729783847], [-1.1442498736307467, 0.36949088558646886], [-1.1441507454679898, 0.3694996552976363], [-1.1440471974014887, 0.369508038985553], [-1.1439392868379046, 0.36951603269357375], [-1.143826962111605, 0.3695236409146822], [-1.143710223846405, 0.3695308637189872], [-1.1435890726877826, 0.36953770117726786]], [0.0, 0.0005538501427314297, 0.0011076681361558054, 0.0016614538131933328, 0.002215207006718231, 0.002768691752159063, 0.003322143707466764, 0.0038755627056006, 0.004428948579475192, 0.004982065253060089 … 104.98257317440857, 104.98265975231972, 104.98275054183942, 104.98284569107507, 104.98294520640296, 104.98304909330444, 104.98315729953856, 104.98326988163824, 104.98338684313389, 104.98350818708334], nothing, DataInterpolations.QuadraticParameterCache{Vector{StaticArraysCore.SVector{2, Float64}}}(StaticArraysCore.SVector{2, Float64}[], StaticArraysCore.SVector{2, Float64}[]), :Forward, DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{Vector{Float64}}([0.0, 0.0005538501427314297, 0.0011076681361558054, 0.0016614538131933328, 0.002215207006718231, 0.002768691752159063, 0.003322143707466764, 0.0038755627056006, 0.004428948579475192, 0.004982065253060089 … 104.98257317440857, 104.98265975231972, 104.98275054183942, 104.98284569107507, 104.98294520640296, 104.98304909330444, 104.98315729953856, 104.98326988163824, 104.98338684313389, 104.98350818708334], Base.RefValue{Int64}(1), true), false, true)
Therefore, we can define the following function to plot the smooth manifold:
using CairoMakie
function manifold_plot(data)
figure = Figure()
axes = Axis(figure[1,1])
for k in eachindex(data)
points = Point2f.(data[k].u)
lines!(axes, points)
end
figure
end
manifold_plot(manifold.data)
Full codes without comments:
using StaticArrays, LinearAlgebra, InvariantManifolds, CairoMakie
function fixedpoints(p)
a , b = p
x1 = (-sqrt(4 * a + b^2 - 2 * b + 1) + b - 1) / (2 * a)
y1 = (1 / 2) * (b^2 / a - b * sqrt(4 * a + b^2 - 2 * b + 1) / a - b / a)
x2 = (sqrt(4 * a + b^2 - 2 * b + 1) + b - 1) / (2 * a)
y2 = (1 / 2) * (b^2 / a + b * sqrt(4 * a + b^2 - 2 * b + 1) / a - b / a)
return SA[x1, y1], SA[x2, y2]
end
function jacobian(x, p)
a, b = p
J = @SMatrix [-2 * a * x[1] 1.0; b 0.0]
return J
end
function henonmap(x, p)
y1 = 1 - p[1] * x[1]^2 + x[2]
y2 = p[2] * x[1]
SA[y1, y2]
end
function henonmap2(x, p)
henonmap(henonmap(x, p), p)
end
para = [1.4, 0.3]
prob = OneDManifoldProblem(henonmap2, para)
saddle = fixedpoints(para)[2]
unstable_direction = eigen(jacobian(fixedpoints([1.4, 0.3])[2], [1.4, 0.3])).vectors[:,1]
segment = gen_segment(saddle, unstable_direction)
manifold = growmanifold(prob, segment, 8)
function manifold_plot(data)
figure = Figure()
axes = Axis(figure[1,1])
for k in eachindex(data)
points = Point2f.(data[k].u)
lines!(axes, points)
end
figure
end
manifold_plot(manifold.data)Oscillator with Periodic Forcing
Now let's consider a higher-order example. Consider the following oscillator with periodic forcing:
\[\begin{aligned} \dot{x}&=y,\\ \dot{y}&=x-\delta x^3+\gamma \cos(\omega t). \end{aligned}\]
When $\gamma=0$, the system has a saddle point at $(0,0)$. After a small periodic perturbation, this saddle point becomes a saddle periodic orbit, which is a saddle point of the mapping $T:X\mapsto \phi(X,2\pi/\omega,0)$, where $\phi(X,t,t_0)$ is the solution of the system under the initial condition $X(t_0)=X\in\mathbb{R}^2$. Fortunately, we can obtain the Jacobian matrix of the mapping $T$ using the solution of the variational equation. The saddle point position and unstable direction of the mapping $T$ can also be obtained through numerical methods.
InvariantManifolds.jl provides a function findsaddle to obtain the saddle point position and unstable direction of $T$. Here's the code demonstrating how to use this function:
using InvariantManifolds, LinearAlgebra, StaticArrays, OrdinaryDiffEq, CairoMakie
f(x, p, t) = SA[x[2], x[1] - p[1]*(x[1]^3) + p[2]*cos(p[3]*t)]
df(x, p, t) = SA[0.0 1.0; 1-p[1]*3*(x[1]^2) 0.0]
initial_guess = SA[0.0, 0.0]
para = [1.0, 0.1, 2.2]
timespan = (0.0, 2pi/para[3])
saddle = findsaddle(f, df, timespan, initial_guess, para)Saddle{2,Float64,Float64}:
saddle: [-0.017123961269067088, 5.508663283951931e-10]
unstable_directions: StaticArraysCore.SVector{2, Float64}[[-0.7071978621024451, -0.7070156885372]]
unstable_eigen_values: [17.380782820480007]The gen_segment function can act directly on the Saddle structure. Therefore, we can use the following code to generate a local manifold:
julia> segment = gen_segment(saddle)50-element Vector{StaticArraysCore.SVector{2, Float64}}: [-0.017123961269067088, 5.508663283951931e-10] [-0.01726828736337371, -0.00014428836516167164] [-0.01741261345768033, -0.00028857728118967163] [-0.017556939551986953, -0.0004328661972176716] [-0.017701265646293575, -0.0005771551132456717] [-0.017845591740600197, -0.0007214440292736718] [-0.01798991783490682, -0.0008657329453016717] [-0.01813424392921344, -0.001010021861329672] [-0.01827857002352006, -0.0011543107773576717] [-0.01842289611782668, -0.0012985996933856717] ⋮ [-0.023041331135638568, -0.005915845006281673] [-0.02318565722994519, -0.006060133922309673] [-0.02332998332425181, -0.006204422838337672] [-0.023474309418558433, -0.006348711754365672] [-0.023618635512865055, -0.006493000670393673] [-0.023762961607171673, -0.006637289586421673] [-0.0239072877014783, -0.006781578502449673] [-0.024051613795784917, -0.006925867418477672] [-0.024195939890091542, -0.0070701563345056725]
Now we can define the nonlinear mapping:
julia> function timeTmap(x, p) prob = ODEProblem{false}(f, x, (0.0, 2pi/p[3]), p) solve(prob, Vern9(), abstol=1e-10)[end] endtimeTmap (generic function with 1 method)
Then create the problem and solve it:
julia> prob = OneDManifoldProblem(timeTmap, para)OneDManifoldProblem: f:Main.timeTmap para: [1.0, 0.1, 2.2] amax: 0.5 d: 0.001 dsmin: 1.0e-6julia> manifold = growmanifold(prob, segment, 7)One-dimensional manifold Curves number: 9 Points number: 361211 Total arc length: 225.98194973455279 Flaw points number: 8 Distance failed points number: 0 Curvature failed points number: 8 Max dα in Flaw Points: 3.937834976139346e-7
Finally, use the function defined in the previous section to plot the results:
manifold_plot(manifold.data)
Full codes without comments:
using StaticArrays, LinearAlgebra, InvariantManifolds, CairoMakie, OrdinaryDiffEq
f(x, p, t) = SA[x[2], x[1] - p[1]*(x[1]^3) + p[2]*cos(p[3]*t)]
df(x, p, t) = SA[0.0 1.0; 1-p[1]*3*(x[1]^2) 0.0]
initial_guess = SA[0.0, 0.0]
para = [1.0, 0.1, 2.2]
timespan = (0.0, 2pi/para[3])
saddle = findsaddle(f, df, timespan, initial_guess, para)
segment = gen_segment(saddle)
function timeTmap(x, p)
prob = ODEProblem{false}(f, x, (0.0, 2pi/p[3]), p)
solve(prob, Vern9(), abstol=1e-10)[end]
end
prob = OneDManifoldProblem(timeTmap, para)
manifold = growmanifold(prob, segment, 7)
function manifold_plot(data)
figure = Figure()
axes = Axis(figure[1,1])
for k in eachindex(data)
points = Point2f.(data[k].u)
lines!(axes, points)
end
figure
end
manifold_plot(manifold.data)