diff --git a/Project.toml b/Project.toml index 35ae220..920d2ae 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PredefinedDynamicalSystems" uuid = "31e2f376-db9e-427a-b76e-a14f56347a14" repo = "https://github.com/JuliaDynamics/PredefinedDynamicalSystems.jl.git" -version = "1.2" +version = "1.3" [deps] DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4" diff --git a/src/continuous_famous_systems.jl b/src/continuous_famous_systems.jl index 46d4507..0b81baf 100644 --- a/src/continuous_famous_systems.jl +++ b/src/continuous_famous_systems.jl @@ -85,7 +85,7 @@ end return SVector{3}(du1, du2, du3) end function chua_jacob(u, p, t) - return SMatrix{3,3}(-p[1]*(1 + chua_element_derivative(u[1], p[3], p[4])),1.0,0.0,p[1],-1.0,-p[2],0.0,1.0,0.0) + return SMatrix{3,3}(-p[1]*(1 + chua_element_derivative(u[1], p[3], p[4])),1.0,0.0,p[1],-1.0,-p[2],0.0,1.0,0.0) end # Helper functions for Chua's circuit. function chua_element(x, m0, m1) @@ -130,7 +130,7 @@ end return SVector{3}(du1, du2, du3) end function roessler_jacob(u, p, t) - return SMatrix{3,3}(0.0,1.0,u[3],-1.0,p[1],0.0,-1.0,0.0,u[1]-p[3]) + return SMatrix{3,3}(0.0,1.0,u[3],-1.0,p[1],0.0,-1.0,0.0,u[1]-p[3]) end """ @@ -170,7 +170,7 @@ end sin_ϕ = sin(φ); cos_ϕ = cos(φ) sin_θ₂ = sin(u[3]); sin_θ₁ = sin(u[1]) Δ = (p[4] + p[5]) - p[5] * cos_ϕ^2 - + du1 = u[2] du2 = (p[5] * (cos_ϕ * (p[2] * u[2]^2 * sin_ϕ + p[1] * sin_θ₂) + p[3] * u[4]^2 * sin_ϕ) - @@ -295,7 +295,7 @@ end `N` is the chain length, `F` the forcing. Jacobian is created automatically. (parameter container only contains `F`) """ -function lorenz96(N::Int, u0 = range(0; length = N, step = 0.1); F=0.01) +function lorenz96(N::Int = 4, u0 = range(0; length = N, step = 0.1); F=0.01) @assert N ≥ 4 "`N` must be at least 4" lor96 = Lorenz96{N}() # create struct return CoupledODEs(lor96, u0, [F]) @@ -398,7 +398,7 @@ function gissinger_rule(u, p, t) return SVector{3}(du1, du2, du3) end function gissinger_jacob(u, p, t) - μ, ν, Γ = p + μ, ν, Γ = p return SMatrix{3,3}(μ,u[3],u[2],-u[3],-ν,u[1],-u[2],u[1],-1) end @@ -433,7 +433,7 @@ function rikitake_jacob(u, p, t) μ, α = p x,y,z = u - return SMatrix{3,3}(-μ,z-α,-y,z,-μ,-x,y,x,0) + return SMatrix{3,3}(-μ,z-α,-y,z,-μ,-x,y,x,0) end """ @@ -472,8 +472,8 @@ function nosehoover_rule(u, p, t) end function nosehoover_jacob(u, p, t) x,y,z = u - - return SMatrix{3,3}(0,-1,0,1,z,-2y,0,y,0) + + return SMatrix{3,3}(0,-1,0,1,z,-2y,0,y,0) end """ @@ -613,7 +613,7 @@ end function ueda_jacob(u, p, t) x,y = u k, B = p - return SMatrix{2,2}(0,-3*x^2,1,-k) + return SMatrix{2,2}(0,-3*x^2,1,-k) end @@ -778,9 +778,9 @@ function stommel_thermohaline_jacob(x, p, t) η1, η2, η3 = p q = abs(T-S) if T ≥ S - return SMatrix{2,2}((-1 - 2T + S), -S,T,(-η3 - T + 2S)) + return SMatrix{2,2}((-1 - 2T + S), -S,T,(-η3 - T + 2S)) else - return SMatrix{2,2}((-1 + 2T - S), S,-T,(-η3 + T - 2S)) + return SMatrix{2,2}((-1 + 2T - S), S,-T,(-η3 + T - 2S)) end end @@ -823,7 +823,7 @@ end end function lorenz84_rule_jacob(u, p, t) F, G, a, b = p - x, y, z = u + x, y, z = u return SMatrix{3,3}(-a,y-b*z,b*y+z,-2y,x-1,b*x,-2z,-b*x,x-1) end @@ -867,7 +867,7 @@ end return SVector{3}(dx, dy, dz) end function lorenzdl_rule_jacob(u, p, t) - x, y, z = u + x, y, z = u return SMatrix{3,3}(-1,-z,y,1,0,x,0,-x,0) end @@ -909,7 +909,7 @@ end """ - kuramoto(D = 20, u0 = range(0, 2π; length = D); + kuramoto(D = 25, u0 = range(0, 2π; length = D); K = 0.3, ω = range(-1, 1; length = D) ) The Kuramoto model[^Kuramoto1975] of `D` coupled oscillators with equation @@ -977,7 +977,7 @@ end function sprott_dissipative_conservative_jacob(u, p, t) a, b, c = p x, y, z = u - + return SMatrix{3,3}(a*y + z,-4x,c - 2x,1 + a*x,b*z,-2y,x,b*y,0) end @@ -1188,7 +1188,7 @@ function hindmarshrose_rule(u, p, t) end function hindmarshrose_jacob(u, p, t) @inbounds begin - a,b,c,d,r,s, xr, I = p + a,b,c,d,r,s, xr, I = p return SMatrix{3,3}(-3*a*u[1]^2 + 2*b*u[1],-2*d*u[1],r*s,1,-1,0,-1,0,-r) end end @@ -1196,7 +1196,7 @@ end """ ```julia hindmarshrose_two_coupled(u0=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6]; -a = 1.0, b = 3.0, d = 5.0, r = 0.001, s = 4.0, xr = -1.6, I = 4.0, +a = 1.0, b = 3.0, c=1.0, d = 5.0, r = 0.001, s = 4.0, xr = -1.6, I = 4.0, k1 = -0.17, k2 = -0.17, k_el = 0.0, xv = 2.0) ``` ```math @@ -1214,7 +1214,7 @@ The default parameter values are taken from article "Dragon-king-like extreme ev coupled bursting neurons", DOI:https://doi.org/10.1103/PhysRevE.97.062311. """ function hindmarshrose_two_coupled(u0=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6]; - a = 1.0, b = 3.0, d = 5.0, r = 0.001, s = 4.0, xr = -1.6, I = 4.0, + a = 1.0, b = 3.0, c=1.0, d = 5.0, r = 0.001, s = 4.0, xr = -1.6, I = 4.0, k1 = -0.17, k2 = -0.17, k_el = 0.0, xv = 2.0) return CoupledODEs(hindmarshrose_coupled_rule, u0, [a, b, c, d, r, s, xr, I, k1, k2, k_el, xv]) end @@ -1225,11 +1225,11 @@ function hindmarshrose_coupled_rule(u, p, t) a, b, c, d, r, s, xr, I, k1, k2, k_el, xv = p x1, y1, z1, x2, y2, z2 = u - du1 = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - vs ) * sigma(x2) + el_link * ( x2 - x1 ) + du1 = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - xv ) * sigma(x2) + k_el * ( x2 - x1 ) du2 = c - d * x1 ^2 - y1 du3 = r * ( s * ( x1 - xr ) - z1 ) - du4 = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - vs ) * sigma(x1) + el_link * ( x1 - x2 ) + du4 = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - xv ) * sigma(x1) + k_el * ( x1 - x2 ) du5 = c - d * x2 ^2 - y2 du6 = r * ( s * ( x2 - xr ) - z2 ) return SVector(du1, du2, du3, du4, du5, du6) @@ -1288,7 +1288,7 @@ end function stuartlandau_jacob(u, p, t) @inbounds begin μ, ω, b = p - + return SMatrix{2,2}(μ - 3*u[1]^2 -u[2]^2 -2*b*u[1]*u[2],-2*u[1]*u[2] +ω +b*u[2]^2 +3*b*u[1]^2, -2*u[1]*u[2] -ω -b*u[1]^2 -3*b*u[2]^2,μ -u[1]^2 -3*u[2]^2 +2*b*u[1]*u[2]) end diff --git a/src/discrete_famous_systems.jl b/src/discrete_famous_systems.jl index dbed941..7b51c40 100644 --- a/src/discrete_famous_systems.jl +++ b/src/discrete_famous_systems.jl @@ -111,7 +111,7 @@ The first `M` entries of the state are the angles, the last `M` are the momenta. """ function coupledstandardmaps end using SparseArrays -function coupledstandardmaps(M::Int, u0 = 0.001rand(2M); +function coupledstandardmaps(M::Int = 2, u0 = 0.001rand(2M); ks = ones(M), Γ = 1.0) SV = SVector{M, Int} @@ -130,7 +130,7 @@ function coupledstandardmaps(M::Int, u0 = 0.001rand(2M); sparseJ = sparse(J) p = vcat(ks, Γ) csm(sparseJ, u0, p, 0) - return DeterministicIteratedMap(csm, u0, p, csm, sparseJ) + return DeterministicIteratedMap(csm, u0, p) end struct CoupledStandardMaps{N} idxs::SVector{N, Int} @@ -148,7 +148,7 @@ function (f::CoupledStandardMaps{N})(xnew::AbstractVector, x, p, n) where {N} xnew[i] = mod2pi(x[i] + xnew[i+N]) end - return xnew + return nothing end function (f::CoupledStandardMaps{M})( J::AbstractMatrix, x, p, n) where {M} @@ -250,23 +250,23 @@ x_n(1 + |2x_n|^{z-1}), & \\quad |x_n| \\le 0.5 \\\\ [2] : Meyer et al., New. J. Phys **20** (2019) """ function pomeau_manneville(u0 = 0.2, z = 2.5) - return DeterministicIteratedMap(pm_rule, u0, [z], pm_jac) + return DeterministicIteratedMap(pm_rule, SVector{1}(u0), [z]) end function pm_rule(x, p, n) - if x < -0.5 - -4x - 3 - elseif -0.5 ≤ x ≤ 0.5 - @inbounds x*(1 + abs(2x)^(p[1]-1)) + if x[1] < -0.5 + SVector{1}(-4x[1] - 3) + elseif -0.5 ≤ x[1] ≤ 0.5 + @inbounds SVector{1}(x[1]*(1 + abs(2x[1])^(p[1]-1))) else - -4x + 3 + SVector{1}(-4x[1] + 3) end end function pomeau_manneville_jacob(x, p, n) - if x < -0.5 + if x[1] < -0.5 -4.0 - elseif -0.5 ≤ x ≤ 0.5 + elseif -0.5 ≤ x[1] ≤ 0.5 @inbounds z = p[1] - 0.5(x^2 * 2^z * (z-1)*abs(x)^(z-3) + 2^z * abs(x)^(z-1) + 2) + 0.5(x[1]^2 * 2^z * (z-1)*abs(x[1])^(z-3) + 2^z * abs(x[1])^(z-1) + 2) else -4.0 end @@ -288,13 +288,13 @@ function's documentation string. [^Manneville1980]: Manneville, P. (1980). Intermittency, self-similarity and 1/f spectrum in dissipative dynamical systems. [Journal de Physique, 41(11), 1235–1243](https://doi.org/10.1051/jphys:0198000410110123500) """ function manneville_simple(x0=0.4; ε = 0.1) - return DeterministicIteratedMap(manneville_f, x0, [ε], manneville_j) + return DeterministicIteratedMap(manneville_f, SVector{1}(x0), [ε]) end function manneville_f(x, p, t) e = p[1] - y = (1+e)*x + (1-e)*x*x - return y%1 + y = (1+e)*x[1] + (1-e)*x[1]*x[1] + return SVector{1}(y%1) end manneville_simple_jacob(x, p, n) = (1+p[1]) + (1-p[1])*2x @@ -411,15 +411,15 @@ function tentmap(u0 = 0.25, μ = 2.0) end function tentmap_rule(x, p, n) μ = p[1] - if x < 0.5 - μ*x + if x[1] < 0.5 + SVector{1}(μ*x[1]) else - μ*(1 - x) + SVector{1}(μ*(1 - x[1])) end end function tentmap_jacob(x, p, n) μ = p[1] - if x < -0.5 + if x[1] < -0.5 μ else -μ @@ -440,14 +440,14 @@ The parameter β controls the dynamics of the map. Its Lyapunov exponent can be At β=2, it becomes the dyadic transformation, also known as the bit shift map, the 2x mod 1 map, the Bernoulli map or the sawtooth map. The typical trajectory for this case is chaotic, though there are countably infinite periodic orbits [^Ott2002]. """ function betatransformationmap(u0 = 0.25; β=2.0) - return DeterministicIteratedMap(betatransformation_rule, u0, [β]) + return DeterministicIteratedMap(betatransformation_rule, SVector{1}(u0), [β]) end function betatransformation_rule(x, p, n) @inbounds β = p[1] - if 0 ≤ x < 1/β - β*x + if 0 ≤ x[1] < 1/β + SVector{1}(β*x[1]) else - β*x - 1 + SVector{1}(β*x - 1) end end function betatransformation_jacob(x, p, n) @@ -520,9 +520,9 @@ end a,b,c,d = p x,y = u t = c - d/(1 + x^2 + y^2) - aux = 2*d/(1+x^2+y^2) - return SMatrix{2,2}(b*(cos(t)-x^2*sin(t)*aux -x*y*cos(t)*aux), b*(sin(t) +x^2*cos(t)*aux)-x*y*sin(t)*aux, - b*(-sin(t) -x*y*sin(t)*aux -y^2*cos(t)*aux), b*(cos(t) -x*y*cos(t)*aux -y^2*sin(t)*aux)) + aux = 2*d/(1+x^2+y^2) + return SMatrix{2,2}(b*(cos(t)-x^2*sin(t)*aux -x*y*cos(t)*aux), b*(sin(t) +x^2*cos(t)*aux)-x*y*sin(t)*aux, + b*(-sin(t) -x*y*sin(t)*aux -y^2*cos(t)*aux), b*(cos(t) -x*y*cos(t)*aux -y^2*sin(t)*aux)) end diff --git a/test/constructors.jl b/test/constructors.jl new file mode 100644 index 0000000..e0b0911 --- /dev/null +++ b/test/constructors.jl @@ -0,0 +1,81 @@ + + +@testset "Discrete Systems Default Constructors" begin + systems = [ + :towel, + :standardmap, + :coupledstandardmaps, + :henon, + :logistic, + :pomeau_manneville, + :manneville_simple, + :arnoldcat, + :grebogi_map, + :nld_coupled_logistic_maps, + :tentmap, + :betatransformationmap, + :rulkovmap, + :ikedamap, + :ulam, + ] + for system in systems + @test @eval PredefinedDynamicalSystems.$system() isa DeterministicIteratedMap + end +end + +@testset "Continous Systems Default Constructors" begin + systems = [ + :lorenz, + :chua, + :roessler, + :double_pendulum, + :henonheiles, + :qbh, + :lorenz96, + :duffing, + :shinriki, + :gissinger, + :rikitake, + :nosehoover, + :antidots, + :ueda, + :magnetic_pendulum, + :fitzhugh_nagumo, + :more_chaos_example, + :thomas_cyclical, + :stommel_thermohaline, + :lorenz84, + :lorenzdl, + :coupled_roessler, + :kuramoto, + :sprott_dissipative_conservative, + :hodgkinhuxley, + :vanderpol, + :lotkavolterra, + :hindmarshrose, + :hindmarshrose_two_coupled, + :stuartlandau_oscillator, + :forced_pendulum, + :riddled_basins, + :morris_lecar, + :sakarya, + :lorenz_bounded, + :swinging_atwood, + :guckenheimer_holmes, + :halvorsen, + :multispecies_competition, + :hyper_roessler, + :hyper_lorenz, + :hyper_qi, + :hyper_jha, + :hyper_wang, + :hyper_xu, + :hyper_bao, + :hyper_cai, + :hyper_lu, + :hyper_pang, + ] + for system in systems + @test @eval PredefinedDynamicalSystems.$system() isa CoupledODEs + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 4c9a8a3..ae9be82 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,3 +2,5 @@ using PredefinedDynamicalSystems using Test @test PredefinedDynamicalSystems.henon() isa DeterministicIteratedMap + +include("constructors.jl")