Merge branch 'master' into spell

NonlinearOscillations · Nov 5, 2024 · e750ace · e750ace
2 parents 4580e71 + ac9278a
commit e750ace
Show file tree

Hide file tree

Showing 43 changed files with 865 additions and 517 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,14 +1,15 @@
 name = "HarmonicBalance"
 uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 authors = ["Jan Kosata <[email protected]>", "Javier del Pino <[email protected]>", "Orjan Ameye <[email protected]>"]
-version = "0.11.1"
+version = "0.12.0"
 
 [deps]
 BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+EndpointRanges = "340492b5-2a47-5f55-813d-aca7ddf97656"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
@@ -42,6 +43,7 @@ DelimitedFiles = "1.9"
 Distances = "0.10.11"
 DocStringExtensions = "0.9.3"
 Documenter = "1.4"
+EndpointRanges = "0.2.2"
 ExplicitImports = "1.6"
 FFTW = "1.8"
 HomotopyContinuation = "2.9"
@@ -51,8 +53,8 @@ Latexify = "0.16"
 ModelingToolkit = "9.34"
 NonlinearSolve = "3.14"
 OrderedCollections = "1.6"
-OrdinaryDiffEqTsit5 = "1.1"
 OrdinaryDiffEqRosenbrock = "1.1"
+OrdinaryDiffEqTsit5 = "1.1"
 Peaks = "0.5"
 Plots = "1.39"
 PrecompileTools = "1.2"
@@ -69,8 +71,8 @@ ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

diff --git a/benchmark/parametron.jl b/benchmark/parametron.jl
@@ -1,4 +1,5 @@
 using HarmonicBalance
+using BenchmarkTools
 
 using Random
 const SEED = 0xd8e5d8df
@@ -18,8 +19,23 @@ dEOM = DifferentialEquation(natural_equation + forces, x)
 add_harmonic!(dEOM, x, ω)
 
 @time harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t);
-@time prob = HarmonicBalance.Problem(harmonic_eq)
 
 fixed = (Ω => 1.0, γ => 1e-2, λ => 5e-2, F => 0, α => 1.0, η => 0.3, θ => 0, ψ => 0)
-varied = ω => range(0.9, 1.1, 20)
-@time res = get_steady_states(prob, varied, fixed; show_progress=false)
+varied = ω => range(0.9, 1.1, 100)
+
+@btime res = get_steady_states(harmonic_eq, WarmUp(), varied, fixed; show_progress=false)
+@btime res = get_steady_states(
+    harmonic_eq, WarmUp(; compile=true), varied, fixed; show_progress=false
+)
+@btime res = get_steady_states(
+    harmonic_eq, Polyhedral(; only_non_zero=true), varied, fixed; show_progress=false
+)
+@btime res = get_steady_states(
+    harmonic_eq, Polyhedral(; only_non_zero=false), varied, fixed; show_progress=false
+)
+@btime res = get_steady_states(
+    harmonic_eq, TotalDegree(; compile=true), varied, fixed; show_progress=false
+)
+@btime res = get_steady_states(
+    harmonic_eq, TotalDegree(), varied, fixed; show_progress=false
+)
diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts
@@ -64,6 +64,7 @@ export default defineConfig({
         text: 'Manual', items: [
           { text: 'Entering equations of motion', link: '/manual/entering_eom.md' },
           { text: 'Computing effective system', link: '/manual/extracting_harmonics' },
+          { text: 'Computing steady states', link: '/manual/methods' },
           { text: 'Krylov-Bogoliubov', link: '/manual/Krylov-Bogoliubov_method' },
           { text: 'Time evolution', link: '/manual/time_dependent' },
           { text: 'Linear response', link: '/manual/linear_response' },
@@ -100,6 +101,7 @@ export default defineConfig({
       "/manual/": [
         { text: 'Entering equations of motion', link: '/manual/entering_eom.md' },
         { text: 'Computing effective system', link: '/manual/extracting_harmonics' },
+        { text: 'Computing steady states', link: '/manual/methods' },
         { text: 'Krylov-Bogoliubov', link: '/manual/Krylov-Bogoliubov_method' },
         { text: 'Time evolution', link: '/manual/time_dependent' },
         { text: 'Linear response', link: '/manual/linear_response' },

diff --git a/docs/src/examples/parametric_via_three_wave_mixing.md b/docs/src/examples/parametric_via_three_wave_mixing.md
@@ -1,5 +1,5 @@
 ```@meta
-EditURL = "parametric_via_three_wave_mixing.jl"
+EditURL = "../../../examples/parametric_via_three_wave_mixing.jl"
 ```
 
 # Parametric Pumping via Three-Wave Mixing
@@ -33,7 +33,7 @@ If we both have quadratic and cubic nonlineariy, we observe the normal duffing o
 varied = (ω => range(0.99, 1.1, 200)) # range of parameter values
 fixed = (α => 1.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025) # fixed parameters
 
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 plot(result; y="u1^2+v1^2")
 ````
 
@@ -43,7 +43,7 @@ If we set the cubic nonlinearity to zero, we recover the driven damped harmonic
 varied = (ω => range(0.99, 1.1, 100))
 fixed = (α => 0.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025)
 
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 plot(result; y="u1^2+v1^2")
 ````
 
@@ -65,17 +65,16 @@ harmonic_eq2 = get_krylov_equations(diff_eq; order=2)
 varied = (ω => range(0.4, 1.1, 500))
 fixed = (α => 1.0, β => 2.0, ω0 => 1.0, γ => 0.001, F => 0.005)
 
-result = get_steady_states(harmonic_eq2, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq2, varied, fixed)
 plot(result; y="v1")
 ````
 
 ````@example parametric_via_three_wave_mixing
 varied = (ω => range(0.4, 0.6, 100), F => range(1e-6, 0.01, 50))
 fixed = (α => 1.0, β => 2.0, ω0 => 1.0, γ => 0.01)
 
-result = get_steady_states(
-    harmonic_eq2, varied, fixed; threading=true, method=:total_degree
-)
+method = TotalDegree()
+result = get_steady_states(harmonic_eq2, method, varied, fixed)
 plot_phase_diagram(result; class="stable")
 ````
 

diff --git a/docs/src/examples/parametron.md b/docs/src/examples/parametron.md
@@ -1,5 +1,5 @@
 ```@meta
-EditURL = "parametron.jl"
+EditURL = "../../../examples/parametron.jl"
 ```
 
 # [Parametrically driven resonator](@id parametron)
@@ -64,7 +64,7 @@ varied = ω => range(0.9, 1.1, 100)
 result = get_steady_states(harmonic_eq, varied, fixed)
 ````
 
-In `get_steady_states`, the default value for the keyword `method=:random_warmup` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `:total_degree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower. The `threading` keyword enables parallel tracking of homotopy paths, and it's set to `false` simply because we are using a single core computer for now.
+In `get_steady_states`, the default method `WarmUp()` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `TotalDegree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower.
 
 After solving the system, we can save the full output of the simulation and the model (e.g. symbolic expressions for the harmonic equations) into a file
 
@@ -100,6 +100,7 @@ The parametrically driven oscillator boasts a stability diagram called "Arnold's
 To perform a 2D sweep over driving frequency $\omega$ and parametric drive strength $\lambda$, we keep `fixed` from before but include 2 variables in `varied`
 
 ````@example parametron
+fixed = (ω₀ => 1.0, γ => 1e-2, F => 1e-3, α => 1.0, η => 0.3)
 varied = (ω => range(0.8, 1.2, 50), λ => range(0.001, 0.6, 50))
 result_2D = get_steady_states(harmonic_eq, varied, fixed);
 nothing #hide

diff --git a/docs/src/examples/wave_mixing.md b/docs/src/examples/wave_mixing.md
@@ -1,5 +1,5 @@
 ```@meta
-EditURL = "wave_mixing.jl"
+EditURL = "../../../examples/wave_mixing.jl"
 ```
 
 # Three Wave Mixing vs four wave mixing
@@ -39,7 +39,7 @@ response with no response at $2\omega$.
 ````@example wave_mixing
 varied = (ω => range(0.9, 1.2, 200)) # range of parameter values
 fixed = (α => 1.0, β => 0.0, ω0 => 1.0, γ => 0.005, F => 0.0025) # fixed parameters
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)# compute steady states
+result = get_steady_states(harmonic_eq, varied, fixed)# compute steady states
 
 p1 = plot(result; y="√(u1^2+v1^2)", legend=:best)
 p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1))
@@ -57,7 +57,7 @@ We would like to investigate the three-wave mixing of the driven Duffing oscilla
 ````@example wave_mixing
 varied = (ω => range(0.9, 1.2, 200))
 fixed = (α => 0.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025)
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 
 p1 = plot(result; y="√(u1^2+v1^2)", legend=:best)
 p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1))
@@ -75,7 +75,7 @@ We would like to investigate the three-wave mixing of the driven Duffing oscilla
 ````@example wave_mixing
 varied = (ω => range(0.9, 1.2, 200))
 fixed = (α => 1.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025)
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 
 p1 = plot(result; y="√(u1^2+v1^2)", legend=:best)
 p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1))

diff --git a/docs/src/manual/methods.md b/docs/src/manual/methods.md
@@ -0,0 +1,18 @@
+# Methods
+
+We offer several methods for solving the nonlinear algebraic equations that arise from the harmonic balance procedure. Each method has different tradeoffs between speed, robustness, and completeness.
+
+## Total Degree Method
+```@docs
+TotalDegree
+```
+
+## Polyhedral Method
+```@docs
+Polyhedral
+```
+
+## Warm Up Method
+```@docs
+WarmUp
+```
diff --git a/docs/src/manual/plotting.md b/docs/src/manual/plotting.md
@@ -12,7 +12,7 @@ The function `plot` is multiple-dispatched to plot 1D and 2D datasets.
 In 1D, the solutions are colour-coded according to the branches obtained by `sort_solutions`. 
 
 ```@docs
-HarmonicBalance.plot(::Result, varags...)
+HarmonicBalance.plot(::HarmonicBalance.Result, varags...)
 ```
 
 

diff --git a/docs/src/manual/solving_harmonics.md b/docs/src/manual/solving_harmonics.md
@@ -8,9 +8,9 @@ having called `get_harmonic_equations`, we need to set all time-derivatives to z
 Once defined, a `Problem` can be solved for a set of input parameters using `get_steady_states` to obtain `Result`.
 
 ```@docs
-Problem
+HarmonicBalance.Problem
 get_steady_states
-Result
+HarmonicBalance.Result
 ```
 
 

diff --git a/docs/src/tutorials/classification.md b/docs/src/tutorials/classification.md
@@ -20,7 +20,7 @@ We performe a 2d sweep in the driving frequency $\omega$ and driving strength $\
 fixed = (ω₀ => 1.0, γ => 0.002, α => 1.0)
 varied = (ω => range(0.99, 1.01, 100), λ => range(1e-6, 0.03, 100))
 
-result_2D = get_steady_states(harmonic_eq, varied, fixed, threading=true)
+result_2D = get_steady_states(harmonic_eq, varied, fixed)
 ```
 By default the steady states of the system are classified by four different catogaries:
 * `physical`: Solutions that are physical, i.e., all variables are purely real.

diff --git a/docs/src/tutorials/limit_cycles.md b/docs/src/tutorials/limit_cycles.md
@@ -62,8 +62,8 @@ Here we reconstruct the results of [Zambon et al., Phys Rev. A 102, 023526 (2020
 ```@example lc
 using HarmonicBalance
 @variables γ F α ω0 F0 η ω J t x(t) y(t);
-eqs = [d(x,t,2) + γ*d(x,t) + ω0^2*x + α*x^3+ 2*J*ω0*(x-y) - F0*cos(ω*t),
-       d(y,t,2) + γ * d(y,t) + ω0^2 * y + α*y^3 + 2*J*ω0*(y-x) - η*F0*cos(ω*t)]
+eqs = [d(x,t,2) + γ*d(x,t) + ω0^2*x + α*x^3 + 2*J*ω0*(x-y) - F0*cos(ω*t),
+       d(y,t,2) + γ*d(y,t) + ω0^2*y + α*y^3 + 2*J*ω0*(y-x) - η*F0*cos(ω*t)]
 diff_eq = DifferentialEquation(eqs, [x,y])
 ```
 
@@ -82,8 +82,7 @@ fixed = (
     J => 154.1e-6,   # coupling term
     α => 3.867e-7,   # Kerr nonlinearity
     ω => 1.4507941,  # pump frequency, resonant with antisymmetric mode (in paper, ħω0 + J)
-    η => -0.08,      # pumping leaking to site 2  (F2 = ηF1)
-    F0 => 0.002       # pump amplitude (overriden in sweeps)
+    η => -0.08     # pumping leaking to site 2  (F2 = ηF1)
 )
 varied = F0 => range(0.002, 0.03, 50)
 

diff --git a/docs/src/tutorials/time_dependent.md b/docs/src/tutorials/time_dependent.md
@@ -65,6 +65,7 @@ plot(time_evo, ["u1", "v1"], harmonic_eq)
 
 Let us compare this to the steady state diagram.
 ```@example time_dependent
+fixed = (ω0 => 1.0, γ => 1e-2, λ => 5e-2, F => 1e-3,  α => 1.0, η => 0.3, θ => 0)
 varied = ω => range(0.9, 1.1, 100)
 result = get_steady_states(harmonic_eq, varied, fixed)
 plot(result, "sqrt(u1^2 + v1^2)")

diff --git a/examples/parametric_via_three_wave_mixing.jl b/examples/parametric_via_three_wave_mixing.jl
@@ -22,15 +22,15 @@ harmonic_eq.equations
 varied = (ω => range(0.99, 1.1, 200)) # range of parameter values
 fixed = (α => 1.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025) # fixed parameters
 
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 plot(result; y="u1^2+v1^2")
 
 # If we set the cubic nonlinearity to zero, we recover the driven damped harmonic oscillator. Indeed, thefirst order the quadratic nonlinearity has no affect on the system.
 
 varied = (ω => range(0.99, 1.1, 100))
 fixed = (α => 0.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025)
 
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 plot(result; y="u1^2+v1^2")
 
 # ## 2nd order Krylov expansion
@@ -50,15 +50,14 @@ harmonic_eq2 = get_krylov_equations(diff_eq; order=2)
 varied = (ω => range(0.4, 1.1, 500))
 fixed = (α => 1.0, β => 2.0, ω0 => 1.0, γ => 0.001, F => 0.005)
 
-result = get_steady_states(harmonic_eq2, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq2, varied, fixed)
 plot(result; y="v1")
 
 #
 
 varied = (ω => range(0.4, 0.6, 100), F => range(1e-6, 0.01, 50))
 fixed = (α => 1.0, β => 2.0, ω0 => 1.0, γ => 0.01)
 
-result = get_steady_states(
-    harmonic_eq2, varied, fixed; threading=true, method=:total_degree
-)
+method = TotalDegree()
+result = get_steady_states(harmonic_eq2, method, varied, fixed)
 plot_phase_diagram(result; class="stable")
diff --git a/examples/parametron.jl b/examples/parametron.jl
@@ -50,7 +50,7 @@ varied = ω => range(0.9, 1.1, 100)
 
 result = get_steady_states(harmonic_eq, varied, fixed)
 
-# In `get_steady_states`, the default value for the keyword `method=:random_warmup` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `:total_degree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower. The `threading` keyword enables parallel tracking of homotopy paths, and it's set to `false` simply because we are using a single core computer for now.
+# In `get_steady_states`, the default method `WarmUp()` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `TotalDegree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower.
 
 # After solving the system, we can save the full output of the simulation and the model (e.g. symbolic expressions for the harmonic equations) into a file
 
@@ -74,6 +74,7 @@ plot(result, "sqrt(u1^2 + v1^2)"; class="all")
 # The parametrically driven oscillator boasts a stability diagram called "Arnold's tongues" delineating zones where the oscillator is stable from those where it is exponentially unstable (if the nonlinearity was absence).  We can retrieve this diagram by calculating the steady states as a function of external detuning $\delta=\omega_L-\omega_0$ and the parametric drive strength $\lambda$.
 
 # To perform a 2D sweep over driving frequency $\omega$ and parametric drive strength $\lambda$, we keep `fixed` from before but include 2 variables in `varied`
+fixed = (ω₀ => 1.0, γ => 1e-2, F => 1e-3, α => 1.0, η => 0.3)
 varied = (ω => range(0.8, 1.2, 50), λ => range(0.001, 0.6, 50))
 result_2D = get_steady_states(harmonic_eq, varied, fixed);
 

diff --git a/examples/wave_mixing.jl b/examples/wave_mixing.jl
@@ -30,7 +30,7 @@ harmonic_eq = get_harmonic_equations(diff_eq)
 
 varied = (ω => range(0.9, 1.2, 200)) # range of parameter values
 fixed = (α => 1.0, β => 0.0, ω0 => 1.0, γ => 0.005, F => 0.0025) # fixed parameters
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)# compute steady states
+result = get_steady_states(harmonic_eq, varied, fixed)# compute steady states
 
 p1 = plot(result; y="√(u1^2+v1^2)", legend=:best)
 p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1))
@@ -46,7 +46,7 @@ plot(p1, p2, p3; layout=(1, 3), size=(900, 300), margin=5mm)
 
 varied = (ω => range(0.9, 1.2, 200))
 fixed = (α => 0.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025)
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 
 p1 = plot(result; y="√(u1^2+v1^2)", legend=:best)
 p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1))
@@ -62,7 +62,7 @@ plot(p1, p2, p3; layout=(1, 3), size=(900, 300), margin=5mm)
 
 varied = (ω => range(0.9, 1.2, 200))
 fixed = (α => 1.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025)
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 
 p1 = plot(result; y="√(u1^2+v1^2)", legend=:best)
 p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1))