@@ -60,8 +60,7 @@ tands = TangentDynamicalSystem(ds; J = towel_jacob)
6060 966.500 μs (10 allocations: 576 bytes) # on my laptop
6161```
6262
63- Here is an example of using [ ` reinit! ` ] ( @ref ) to efficiently iterate over different parameter values, and parallelize via ` Threads ` , to compute the exponents over a given parameter range.
64-
63+ Here is an example of using [ ` reinit! ` ] ( @ref ) to iterate over different parameter values to compute the exponents over a given parameter range.
6564
6665``` @example MAIN
6766using ChaosTools, CairoMakie
@@ -74,14 +73,9 @@ tands = TangentDynamicalSystem(ds; J = henon_jacob)
7473as = 0.8:0.005:1.225;
7574λs = zeros(length(as), 2)
7675
77- # Since `DynamicalSystem`s are mutable, we need to copy to parallelize
78- systems = [deepcopy(tands) for _ in 1:Threads.nthreads()-1]
79- pushfirst!(systems, tands)
80-
81- Threads.@threads for i in eachindex(as)
82- system = systems[Threads.threadid()]
83- set_parameter!(system, 1, as[i])
84- λs[i, :] .= lyapunovspectrum(system, 10000; Ttr = 500)
76+ for i in eachindex(as)
77+ set_parameter!(ds, 1, as[i])
78+ λs[i, :] .= lyapunovspectrum(ds, 10000; Ttr = 500)
8579end
8680
8781fig = Figure()
9286fig
9387```
9488
89+ You can parallelize this code as well, see the main tutorial of DynamicalSystems.jl for an example.
90+
9591## Maximum Lyapunov Exponent
9692
9793It is possible to get only the maximum Lyapunov exponent simply by giving
@@ -216,7 +212,7 @@ pidx = 3 # set to dummy, not used anywhere (no drift)
216212ρ,times = ensemble_averaged_pairwise_distance(ds,init_states,100,pidx;Ttr=5000)
217213λ_inst = lyapunov_instant(ρ,times;interval=20:30) #fit to middle part of the curve (slope is constant until saturation)
218214λ = lyapunov(ds,1000;Ttr=5000) #standard (Benettin) way
219- @show λ_inst, λ
215+ @show λ_inst, λ
220216```
221217
222218Now look at the nonautonomous Duffing map with drifting ε parameter:
@@ -236,17 +232,17 @@ end
236232 return SVector(dx1, dx2)
237233end
238234
239- duffing = duffing_drift()
235+ duffing = duffing_drift()
240236duffing_map = StroboscopicMap(duffing,2π)
241- init_states = randn(5000,2) #use an ensemble of 5000
237+ init_states = randn(5000,2) #use an ensemble of 5000
242238pidx = 4 #ε is the fourth parameter
243239ρ,times = ensemble_averaged_pairwise_distance(duffing_map,StateSpaceSet(init_states),100,pidx;Ttr=20)
244240
245- #measure slope of ρ at two places
241+ #measure slope of ρ at two places
246242λ_inst = lyapunov_instant(ρ,times;interval=5:10)
247243λ_inst2 = lyapunov_instant(ρ,times;interval=22:25)
248244
249- fig,ax,obj = scatter(times, ρ,
245+ fig,ax,obj = scatter(times, ρ,
250246 markersize = 6,
251247 color = :gray10,
252248 label = L"\omega = 1, \beta = 0.2, \epsilon_0 = 0.4, \alpha=0.00045",
@@ -258,7 +254,7 @@ lines!(ax, times[22:25], ρ[22] .+ λ_inst2*[0:3;],color = (:red, 0.7), linewidt
258254text!(ax, times[5]+8, ρ[10],text = L"\lambda = %$(round(λ_inst;digits=3))",
259255color = :red,align = (:left, :center),fontsize = 18)
260256
261- text!(ax, times[22]+8, ρ[24],text = L"\lambda = %$(round(λ_inst2;digits=3))",
257+ text!(ax, times[22]+8, ρ[24],text = L"\lambda = %$(round(λ_inst2;digits=3))",
262258color = :red, align = (:left, :center), fontsize = 18)
263259
264260axislegend(ax, position = :rb, nbanks = 2,labelsize = 18)
0 commit comments