Merge pull request #79 from anandijain/aj/attractors

add common strange attractor systems and problems

Author: ChrisRackauckas

Project.toml

@@ -1,22 +1,23 @@
 name = "DiffEqProblemLibrary"
 uuid = "a077e3f3-b75c-5d7f-a0c6-6bc4c8ec64a9"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "4.13.1"
+version = "4.13.2"
 
 [deps]
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqOperators = "9fdde737-9c7f-55bf-ade8-46b3f136cc48"
+Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
-Catalyst = "5.1, 6.0"
+Catalyst = "9"
 DiffEqBase = "6"
 DiffEqOperators = "4"
-ModelingToolkit = "4, 5.0"
+ModelingToolkit = "6"
 julia = "1"
 
 [extras]

src/DiffEqProblemLibrary.jl

@@ -1,6 +1,7 @@
 module DiffEqProblemLibrary
 
 module ODEProblemLibrary
+using Latexify
 function importodeproblems()
   @isdefined(prob_ode_linear) ||
   include(joinpath(@__DIR__, "ode/ode_premade_problems.jl"));

src/jump_premade_problems.jl

@@ -137,7 +137,7 @@ end kon kAon koff kAoff kAp kAdp
 rsi   = rates_sym_to_idx
 rates = params[[rsi[:kon], rsi[:kAon], rsi[:koff], rsi[:kAoff], rsi[:kAp], rsi[:kAdp]]]
 u0    = zeros(Int,9)
-statesyms = ModelingToolkit.tosymbol.(ModelingToolkit.operation.(rs.states))
+statesyms = ModelingToolkit.tosymbol.(ModelingToolkit.operation.(states(rs)))
 u0[ findfirst(isequal(:S1), statesyms)] = params[1]
 u0[ findfirst(isequal(:S2), statesyms)] = params[2]
 u0[ findfirst(isequal(:S3), statesyms)] = params[3]
@@ -194,10 +194,10 @@ function construct_genenetwork(N)
     genenetwork
 end
 rs = construct_genenetwork(N)
-u0 = zeros(Int, length(rs.states))
-statesyms = ModelingToolkit.tosymbol.(ModelingToolkit.operation.(rs.states))
+u0 = zeros(Int, length(states(rs)))
+statesyms = ModelingToolkit.tosymbol.(ModelingToolkit.operation.(states(rs)))
 for i = 1:(2*N)
-    u0[findfirst(isequal(G[i]),rs.states)] = 1
+    u0[findfirst(isequal(G[i]),states(rs))] = 1
 end
 tf = 2000.0
 prob = DiscreteProblem(rs, u0, (0.0, tf))

src/ode/ode_premade_problems.jl

@@ -5,17 +5,22 @@ Random.seed!(100)
 
 #ODE Example Problems
 export prob_ode_linear, prob_ode_bigfloatlinear, prob_ode_2Dlinear,
-       prob_ode_large2Dlinear, prob_ode_bigfloat2Dlinear, prob_ode_rigidbody,
-       prob_ode_2Dlinear_notinplace, prob_ode_vanderpol, prob_ode_vanderpol_stiff,
-       prob_ode_lorenz, prob_ode_rober, prob_ode_threebody, prob_ode_mm_linear, prob_ode_pleiades,
-       prob_ode_brusselator_1d, prob_ode_brusselator_2d, prob_ode_orego,
-       prob_ode_hires, prob_ode_pollution, prob_ode_filament,
-       prob_ode_nonlinchem,
-       SolverDiffEq, filament_prob
+        prob_ode_large2Dlinear, prob_ode_bigfloat2Dlinear, prob_ode_rigidbody,
+        prob_ode_2Dlinear_notinplace, prob_ode_vanderpol, prob_ode_vanderpol_stiff,
+        prob_ode_rober, prob_ode_threebody, prob_ode_mm_linear, prob_ode_pleiades,
+        prob_ode_brusselator_1d, prob_ode_brusselator_2d, prob_ode_orego,
+        prob_ode_hires, prob_ode_pollution, prob_ode_filament,
+        prob_ode_nonlinchem,
+        SolverDiffEq, filament_prob,
+        thomas, lorenz, aizawa, dadras, chen, rossler, rabinovich_fabrikant, sprott, hindmarsh_rose,
+        prob_ode_thomas, prob_ode_lorenz, prob_ode_aizawa, prob_ode_dadras, 
+        prob_ode_chen, prob_ode_rossler, prob_ode_rabinovich_fabrikant, prob_ode_sprott,
+        prob_ode_hindmarsh_rose
 
 include("ode_linear_prob.jl")
 include("ode_simple_nonlinear_prob.jl")
 include("brusselator_prob.jl")
 include("pollution_prob.jl")
 include("filament_prob.jl")
 include("nonlinchem.jl")
+include("strange_attractors.jl")

src/ode/ode_simple_nonlinear_prob.jl

@@ -51,10 +51,10 @@ prob_ode_fitzhughnagumo = ODEProblem(fitz,[1.0;1.0],(0.0,1.0),(0.7,0.8,1/12.5,0.
 #Van der Pol Equations
 @parameters t μ
 @variables x(t) y(t)
-@derivatives D'~t
+D = Differential(t)
 eqs = [D(y) ~ μ*((1-x^2)*y - x),
        D(x) ~ y]
-de = ODESystem(eqs)
+de = ODESystem(eqs; name=:van_der_pol)
 van = ODEFunction(de, [y,x], [μ], jac=true)
 
 """
@@ -87,16 +87,16 @@ with ``μ=10^6`` and ``u_0=[0,\\sqrt{3}]``
 
 Stiff parameters.
 """
-prob_ode_vanstiff = ODEProblem(van,[0;sqrt(3)],(0.0,1.0),1e6)
+prob_ode_vanderpol_stiff = ODEProblem(van,[0;sqrt(3)],(0.0,1.0),1e6)
 
 # ROBER
 @parameters t k₁ k₂ k₃
 @variables y₁(t) y₂(t) y₃(t)
-@derivatives D'~t
+D = Differential(t)
 eqs = [D(y₁) ~ -k₁*y₁+k₃*y₂*y₃,
        D(y₂) ~ k₁*y₁-k₂*y₂^2-k₃*y₂*y₃,
        D(y₃) ~ k₂*y₂^2]
-de = ODESystem(eqs)
+de = ODESystem(eqs; name=:rober)
 rober = ODEFunction(de, [y₁,y₂,y₃], [k₁,k₂,k₃], jac=true)
 
 """
@@ -169,11 +169,11 @@ prob_ode_threebody = ODEProblem(threebody,[0.994, 0.0, 0.0, big(-2.0015851063790
 
 @parameters t I₁ I₂ I₃
 @variables y₁(t) y₂(t) y₃(t)
-@derivatives D'~t
+D = Differential(t)
 eqs = [D(y₁) ~ I₁*y₂*y₃,
        D(y₂) ~ I₂*y₁*y₃,
        D(y₃) ~ I₃*y₁*y₂]
-de = ODESystem(eqs)
+de = ODESystem(eqs; name=:rigid_body)
 rigid = ODEFunction(de, [y₁,y₂,y₃], [I₁,I₂,I₃], jac=true)
 
 """
@@ -321,7 +321,7 @@ prob_ode_mm_linear = ODEProblem(mm_f,rand(4),(0.0,1.0))
 
 @parameters t p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12
 @variables y1(t) y2(t) y3(t) y4(t) y5(t) y6(t) y7(t) y8(t)
-@derivatives D'~t
+D = Differential(t)
 eqs = [D(y1) ~ -p1*y1 + p2*y2 + p3*y3 + p4,
        D(y2) ~ p1*y1 - p5*y2,
        D(y3) ~ -p6*y3 + p2*y4 + p7*y5,
@@ -331,7 +331,7 @@ eqs = [D(y1) ~ -p1*y1 + p2*y2 + p3*y3 + p4,
                 p2*y6 + p11*y7,
        D(y7) ~  p10*y6*y8 - p12*y7,
        D(y8) ~ -p10*y6*y8 + p12*y7]
-de = ODESystem(eqs)
+de = ODESystem(eqs; name=:hires)
 hires = ODEFunction(de, [y1,y2,y3,y4,y5,y6,y7,y8],
                         [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12],
                         jac=true)
@@ -368,11 +368,11 @@ prob_ode_hires = ODEProblem(hires,u0,(0.0,321.8122), (1.71, 0.43, 8.32, 0.0007,
 
 @parameters t p1 p2 p3
 @variables y1(t) y2(t) y3(t)
-@derivatives D'~t
+D = Differential(t)
 eqs = [D(y1) ~ p1*(y2+y1*(1-p2*y1-y2)),
        D(y2) ~ (y3-(1+y1)*y2)/p1,
        D(y3) ~ p3*(y1-y3)]
-de = ODESystem(eqs)
+de = ODESystem(eqs; name=:orego)
 jac = calculate_jacobian(de)
 orego = ODEFunction(de, [y1,y2,y3], [p1,p2,p3], jac=true)
 

src/ode/strange_attractors.jl

@@ -0,0 +1,221 @@
+# Strange attractors from https://www.dynamicmath.xyz/strange-attractors/ and
+# https://en.wikipedia.org/wiki/List_of_chaotic_maps
+# Opted for the equations as reported in papers 
+
+# Thomas
+@parameters t b=0.208186
+@variables x(t)=1 y(t)=0 z(t)=0
+D = Differential(t)
+
+eqs = [D(x) ~ sin(y) - b*x,
+       D(y) ~ sin(z) - b*y,
+       D(z) ~ sin(x) - b*z]
+
+@named thomas = ODESystem(eqs)
+
+"""
+Thomas' cyclically symmetric attractor equations
+
+```math
+$(latexify(thomas))
+```
+
+[Reference](https://www.worldscientific.com/doi/abs/10.1142/S0218127499001383)
+
+[Wikipedia](https://en.wikipedia.org/wiki/Thomas%27_cyclically_symmetric_attractor)
+"""
+prob_ode_thomas= ODEProblem(thomas,[],(0.0,1.0))
+
+
+# Lorenz
+@parameters t σ=10 ρ=28 β=8/3
+@variables x(t)=1 y(t)=0 z(t)=0
+D = Differential(t)
+
+eqs = [D(x) ~ σ*(y-x),
+       D(y) ~ x*(ρ-z)-y,
+       D(z) ~ x*y - β*z]
+
+@named lorenz = ODESystem(eqs)
+
+"""
+Lorenz equations
+
+```math
+$(latexify(lorenz))
+```
+
+[Reference](https://journals.ametsoc.org/view/journals/atsc/20/2/1520-0469_1963_020_0130_dnf_2_0_co_2.xml)
+
+[Wikipedia](https://en.wikipedia.org/wiki/Lorenz_system)
+"""
+prob_ode_lorenz = ODEProblem(lorenz,[],(0.0,1.0))
+
+
+# Aizawa
+@parameters t a=0.95 b=0.7 c=0.6 d=3.5 e=0.25 f=0.1
+@variables x(t)=1 y(t)=0 z(t)=0
+D = Differential(t)
+
+eqs = [D(x) ~ (z-b)*x-d*y,
+       D(y) ~ d*x+(z-b)*y,
+       D(z) ~ c+a*z-z^3/3-(x^2+y^2)*(1+e*z)+ f*z*x^3]
+
+@named aizawa = ODESystem(eqs)
+
+"""
+Aizawa equations
+
+```math
+$(latexify(aizawa))
+```
+
+[Reference](https://journals.ametsoc.org/view/journals/atsc/20/2/1520-0469_1963_020_0130_dnf_2_0_co_2.xml)
+
+"""
+prob_ode_aizawa = ODEProblem(aizawa,[],(0.0,1.0))
+
+
+# Dadras
+@parameters t a=3 b=2.7 c=1.7 d=2 e=9
+@variables x(t)=1 y(t)=0 z(t)=0
+D = Differential(t)
+
+eqs = [D(x) ~ y-a*x+b*y*z,
+       D(y) ~ c*y-x*z+z,
+       D(z) ~ d*x*y-e*z]
+
+@named dadras = ODESystem(eqs)
+
+"""
+Dadras equations
+
+```math
+$(latexify(dadras))
+```
+
+[Reference](https://www.sciencedirect.com/science/article/abs/pii/S0375960109009591)
+
+"""
+prob_ode_dadras = ODEProblem(dadras,[],(0.0,1.0))
+
+
+# chen
+@parameters t a=35 b=3 c=28
+@variables x(t)=1 y(t)=0 z(t)=0
+D = Differential(t)
+
+eqs = [D(x) ~ a*(y-x),
+       D(y) ~ (c-a)*x-x*z+c*y,
+       D(z) ~ x*y-b*z]
+
+@named chen = ODESystem(eqs)
+
+"""
+chen equations
+
+```math
+$(latexify(chen))
+```
+
+[Reference](https://www.worldscientific.com/doi/abs/10.1142/S0218127499001024)
+
+"""
+prob_ode_chen = ODEProblem(chen,[],(0.0,1.0))
+
+
+# rossler
+@parameters t a=0.2 b=0.2 c=5.7
+@variables x(t)=1 y(t)=0 z(t)=0
+D = Differential(t)
+
+eqs = [D(x) ~ -(y+z),
+       D(y) ~ x+a*y,
+       D(z) ~ b+z*(x-c)]
+
+@named rossler = ODESystem(eqs)
+
+"""
+rossler equations
+
+```math
+$(latexify(rossler))
+```
+
+[Reference](https://www.sciencedirect.com/science/article/abs/pii/0375960176901018)
+[Wikipedia](https://en.wikipedia.org/wiki/R%C3%B6ssler_attractor)
+
+"""
+prob_ode_rossler = ODEProblem(rossler,[],(0.0,1.0))
+
+
+# rabinovich_fabrikant
+@parameters t a=0.14 b=0.10
+@variables x(t)=1 y(t)=0 z(t)=0
+D = Differential(t)
+
+eqs = [D(x) ~ y*(z-1+x^2) + b*x,
+       D(y) ~ x*(3*z+1-x^2)+b*y,
+       D(z) ~ -2*z*(a+x*y)]
+
+@named rabinovich_fabrikant = ODESystem(eqs)
+
+"""
+rabinovich_fabrikant equations
+
+```math
+$(latexify(rabinovich_fabrikant))
+```
+
+[Reference](https://en.wikipedia.org/wiki/Rabinovich%E2%80%93Fabrikant_equations)
+
+"""
+prob_ode_rabinovich_fabrikant = ODEProblem(rabinovich_fabrikant,[],(0.0,1.0))
+
+
+# sprott
+@parameters t a=2.07 b=1.79
+@variables x(t)=1 y(t)=0 z(t)=0
+D = Differential(t)
+
+eqs = [D(x) ~ y+a*x*y+x*z,
+       D(y) ~ 1-b*x^2+y*z,
+       D(z) ~ x-x^2-y^2]
+
+@named sprott = ODESystem(eqs)
+
+"""
+sprott equations
+
+```math
+$(latexify(sprott))
+```
+
+[Reference](http://sprott.physics.wisc.edu/pubs/paper423.pdf)
+
+"""
+prob_ode_sprott = ODEProblem(sprott,[],(0.0,1.0))
+
+
+# hindmarsh_rose
+@parameters t a=1 b=3 c=1 d=5 r=1e-2 s=4 xr=-8/5 i=5
+@variables x(t)=1 y(t)=0 z(t)=0
+D = Differential(t)
+
+eqs = [D(x) ~ y-a*x^3+b*x^2-z+i,
+       D(y) ~ c-d*x^2-y,
+       D(z) ~ r*(s*(x-xr)-z)]
+
+@named hindmarsh_rose = ODESystem(eqs)
+
+"""
+hindmarsh_rose equations
+
+```math
+$(latexify(hindmarsh_rose))
+```
+
+[Reference](https://en.wikipedia.org/wiki/Hindmarsh%E2%80%93Rose_model)
+
+"""
+prob_ode_hindmarsh_rose = ODEProblem(hindmarsh_rose,[],(0.0,1.0))