diff --git a/Project.toml b/Project.toml index ae458653..57273395 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ 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 index 4b257a8c..d1d23ea3 100644 --- a/benchmark/parametron.jl +++ b/benchmark/parametron.jl @@ -1,4 +1,5 @@ using HarmonicBalance +using BenchmarkTools using Random const SEED = 0xd8e5d8df @@ -18,8 +19,13 @@ 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 index b823cdcc..6fc8ed59 100644 --- 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 index 7908825c..3bc2b036 100644 --- 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,7 +65,7 @@ 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") ```` @@ -73,9 +73,8 @@ 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/docs/src/examples/parametron.md b/docs/src/examples/parametron.md index 26b5c929..4db5b685 100644 --- 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 index 95df2bad..2be86998 100644 --- 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 new file mode 100644 index 00000000..80e8b52a --- /dev/null +++ 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 +``` \ No newline at end of file diff --git a/docs/src/manual/plotting.md b/docs/src/manual/plotting.md index bb096aab..6db443aa 100644 --- 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 index 607d878a..65057e61 100644 --- 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 index 8996b79a..b79178ec 100644 --- 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 index bdeaab84..acd790a9 100644 --- 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 index 6065960a..0bd13c9f 100644 --- 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 index 3b86887c..4a084fea 100644 --- a/examples/parametric_via_three_wave_mixing.jl +++ b/examples/parametric_via_three_wave_mixing.jl @@ -22,7 +22,7 @@ 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. @@ -30,7 +30,7 @@ plot(result; y="u1^2+v1^2") 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,7 +50,7 @@ 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") # @@ -58,7 +58,6 @@ 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 index e9cff811..0db5a953 100644 --- 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 index 2af140f2..db548ac0 100644 --- 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)) diff --git a/src/DifferentialEquation.jl b/src/DifferentialEquation.jl index 01235ddc..01721e9e 100644 --- a/src/DifferentialEquation.jl +++ b/src/DifferentialEquation.jl @@ -1,3 +1,81 @@ +""" +$(TYPEDEF) + +Holds differential equation(s) of motion and a set of harmonics to expand each variable. +This is the primary input for `HarmonicBalance.jl` ; after inputting the equations, the harmonics + ansatz needs to be specified using `add_harmonic!`. + +# Fields +$(TYPEDFIELDS) + +## Example +```julia-repl +julia> @variables t, x(t), y(t), ω0, ω, F, k; + +# equivalent ways to enter the simple harmonic oscillator +julia> DifferentialEquation(d(x,t,2) + ω0^2 * x - F * cos(ω*t), x); +julia> DifferentialEquation(d(x,t,2) + ω0^2 * x ~ F * cos(ω*t), x); + +# two coupled oscillators, one of them driven +julia> DifferentialEquation([d(x,t,2) + ω0^2 * x - k*y, d(y,t,2) + ω0^2 * y - k*x] .~ [F * cos(ω*t), 0], [x,y]); +``` +""" +mutable struct DifferentialEquation + """Assigns to each variable an equation of motion.""" + equations::OrderedDict{Num,Equation} + """Assigns to each variable a set of harmonics.""" + harmonics::OrderedDict{Num,OrderedSet{Num}} + + function DifferentialEquation(eqs) + return new(eqs, OrderedDict(var => OrderedSet() for var in keys(eqs))) + end + + # uses the above constructor if no harmonics defined + function DifferentialEquation(eqs::Vector{Equation}, vars::Vector{Num}) + return DifferentialEquation(OrderedDict(zip(vars, eqs))) + end + + # if expressions are entered instead of equations, automatically set them = 0 + function DifferentialEquation(exprs::Vector{Num}, vars::Vector{Num}) + return DifferentialEquation(exprs .~ Int(0), vars) + end + + function DifferentialEquation(eq::Equation, var::Num) + typerhs = typeof(eq.rhs) + typelhs = typeof(eq.lhs) + if eq.rhs isa AbstractVector || eq.lhs isa AbstractVector + throw( + ArgumentError( + "The equation is of the form $(typerhs)~$(typelhs) is not supported. Commenly one forgot to broadcast the equation symbol `~`.", + ), + ) + end + return DifferentialEquation([eq], [var]) + end + function DifferentialEquation(eq::Equation, vars::Vector{Num}) + typerhs = typeof(eq.rhs) + typelhs = typeof(eq.lhs) + throw( + ArgumentError( + "The variables are of type $(typeof(vars)). Whereas you gave one equation of type $(typerhs)~$(typelhs). Commenly one forgot to broadcast the equation symbol `~`.", + ), + ) + end + DifferentialEquation(lhs::Num, var::Num) = DifferentialEquation([lhs ~ Int(0)], [var]) +end + +function Base.show(io::IO, diff_eq::DifferentialEquation) + println(io, "System of ", length(keys(diff_eq.equations)), " differential equations") + println(io, "Variables: ", join(keys(diff_eq.equations), ", ")) + print(io, "Harmonic ansatz: ") + for var in keys(diff_eq.harmonics) + print(io, string(var), " => ", join(string.(diff_eq.harmonics[var]), ", ")) + print(io, "; ") + end + println(io, "\n") + return [println(io, eq) for eq in values(diff_eq.equations)] +end + """ $(TYPEDSIGNATURES) Add the harmonic `ω` to the harmonic ansatz used to expand the variable `var` in `diff_eom`. diff --git a/src/HarmonicBalance.jl b/src/HarmonicBalance.jl index 9d0cb4fd..f8d22f61 100644 --- a/src/HarmonicBalance.jl +++ b/src/HarmonicBalance.jl @@ -13,6 +13,7 @@ using OrderedCollections: OrderedDict, OrderedSet using ProgressMeter: ProgressMeter, Progress using LinearAlgebra: eigvals using Random: Random # for setting seed +using EndpointRanges: EndpointRanges using Distances: Distances using BijectiveHilbert: BijectiveHilbert, Simple2D, decode_hilbert!, encode_hilbert @@ -40,21 +41,25 @@ using .ExprUtils: is_harmonic, substitute_all, drop_powers, count_derivatives include("extention_functions.jl") include("utils.jl") include("types.jl") - include("DifferentialEquation.jl") include("HarmonicVariable.jl") include("HarmonicEquation.jl") +include("Problem.jl") +include("Result.jl") +include("methods.jl") + include("solve_homotopy.jl") include("sorting.jl") include("classification.jl") + include("saving.jl") include("transform_solutions.jl") include("plotting_Plots.jl") -export show, *, @variables, d, ComplexF64, Float64, IM_TOL +export show, *, @variables, d, IM_TOL +export WarmUp, TotalDegree, Polyhedral -export ParameterRange, ParameterList, StateDict, SteadyState, ParameterVector -export DifferentialEquation, HarmonicVariable, HarmonicEquation, Problem, Result +export DifferentialEquation, HarmonicVariable, HarmonicEquation export get_steady_states, get_single_solution, get_harmonic_equations, add_harmonic! export get_variables, get_independent_variables, classify_branch, classify_solutions! export rearrange_standard @@ -86,14 +91,14 @@ using .FFTWExt using PrecompileTools: @setup_workload, @compile_workload -@setup_workload begin - # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the - # precompile file and potentially make loading faster. - @compile_workload begin - # all calls in this block will be precompiled, regardless of whether - # they belong to your package or not (on Julia 1.8 and higher) - include("precompilation.jl") - end -end +# @setup_workload begin +# # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the +# # precompile file and potentially make loading faster. +# @compile_workload begin +# # all calls in this block will be precompiled, regardless of whether +# # they belong to your package or not (on Julia 1.8 and higher) +# include("precompilation.jl") +# end +# end end # module diff --git a/src/HarmonicEquation.jl b/src/HarmonicEquation.jl index 7822e915..1a126851 100644 --- a/src/HarmonicEquation.jl +++ b/src/HarmonicEquation.jl @@ -1,3 +1,57 @@ + +""" +$(TYPEDEF) + +Holds a set of algebraic equations governing the harmonics of a `DifferentialEquation`. + +# Fields +$(TYPEDFIELDS) +""" +mutable struct HarmonicEquation + """A set of equations governing the harmonics.""" + equations::Vector{Equation} + """A set of variables describing the harmonics.""" + variables::Vector{HarmonicVariable} + """The parameters of the equation set.""" + parameters::Vector{Num} + "The natural equation (before the harmonic ansatz was used)." + natural_equation::DifferentialEquation + + # use a self-referential constructor with _parameters + function HarmonicEquation(equations, variables, nat_eq) + return (x = new(equations, variables, Vector{Num}([]), nat_eq); + x.parameters = _parameters(x); + x) + end + function HarmonicEquation(equations, variables, parameters, natural_equation) + return new(equations, variables, parameters, natural_equation) + end +end + +function Base.show(io::IO, eom::HarmonicEquation) + println(io, "A set of ", length(eom.equations), " harmonic equations") + println(io, "Variables: ", join(string.(get_variables(eom)), ", ")) + println(io, "Parameters: ", join(string.(eom.parameters), ", ")) + println(io, "\nHarmonic ansatz: ", _show_ansatz(eom)) + println(io, "\nHarmonic equations:") + return [println(io, "\n", eq) for eq in eom.equations] +end + +"""Gives the full harmonic ansatz used to construct `eom`.""" +function _show_ansatz(eom::HarmonicEquation) + output = "" + vars = unique(getfield.(eom.variables, :natural_variable)) + for nat_var in vars + # the Hopf variable (limit cycle frequency) does not contribute a term + harm_vars = filter( + x -> isequal(nat_var, x.natural_variable) && x.type !== "Hopf", eom.variables + ) + ansatz = join([_show_ansatz(var) for var in harm_vars], " + ") + output *= "\n" * string(nat_var) * " = " * ansatz + end + return output +end + Base.show(eom::HarmonicEquation) = show_fields(eom) """ @@ -140,9 +194,6 @@ function Symbolics.get_variables(eom::HarmonicEquation)::Vector{Num} return get_variables.(eom.variables) end -Symbolics.get_variables(p::Problem)::Vector{Num} = get_variables(p.eom) -Symbolics.get_variables(res::Result)::Vector{Num} = get_variables(res.problem) - "Get the parameters (not time nor variables) of a HarmonicEquation" function _parameters(eom::HarmonicEquation) all_symbols = flatten([ diff --git a/src/HarmonicVariable.jl b/src/HarmonicVariable.jl index 332f987a..ebed7618 100644 --- a/src/HarmonicVariable.jl +++ b/src/HarmonicVariable.jl @@ -1,3 +1,44 @@ +""" +$(TYPEDEF) + +Holds a variable stored under `symbol` describing the harmonic `ω` of `natural_variable`. + +# Fields +$(TYPEDFIELDS) +""" +mutable struct HarmonicVariable + """Symbol of the variable in the HarmonicBalance namespace.""" + symbol::Num + """Human-readable labels of the variable, used for plotting.""" + name::String + """Type of the variable (u or v for quadratures, a for a constant, Hopf for Hopf etc.)""" + type::String + """The harmonic being described.""" + ω::Num + """The natural variable whose harmonic is being described.""" + natural_variable::Num +end + +function Base.show(io::IO, hv::HarmonicVariable) + return println( + io, + "Harmonic variable ", + string.(hv.symbol) * " for harmonic ", + string(hv.ω), + " of ", + string(hv.natural_variable), + ) +end + +"""Gives the relation between `var` and the underlying natural variable.""" +function _show_ansatz(var::HarmonicVariable) + t = var.natural_variable.val.arguments + t = length(t) == 1 ? string(t[1]) : error("more than 1 independent variable") + ω = string(var.ω) + terms = Dict("u" => "*cos(" * ω * t * ")", "v" => "*sin(" * ω * t * ")", "a" => "") + return string(string(var.symbol) * terms[var.type]) +end + # pretty-printing Base.display(var::HarmonicVariable) = display(var.name) Base.display(var::Vector{HarmonicVariable}) = display.(getfield.(var, Symbol("name"))) diff --git a/src/LimitCycles/LimitCycles.jl b/src/LimitCycles/LimitCycles.jl index 1babaf59..8e6ba97a 100644 --- a/src/LimitCycles/LimitCycles.jl +++ b/src/LimitCycles/LimitCycles.jl @@ -6,6 +6,11 @@ using Symbolics: Symbolics, Num, expand_derivatives using HarmonicBalance using HarmonicBalance: + HarmonicBalanceMethod, + WarmUp, + Problem, + Result, + get_steady_states, order_branches!, find_branch_order, _remove_brackets, @@ -15,6 +20,7 @@ using HarmonicBalance: _is_physical, substitute_all, var_name + using HarmonicBalance.LinearResponse: get_implicit_Jacobian using HarmonicBalance.ExprUtils: get_all_terms diff --git a/src/LimitCycles/gauge_fixing.jl b/src/LimitCycles/gauge_fixing.jl index f84f3936..b579c742 100644 --- a/src/LimitCycles/gauge_fixing.jl +++ b/src/LimitCycles/gauge_fixing.jl @@ -85,16 +85,7 @@ function _choose_fixed(eom, ω_lc) return first(vars) # This is arbitrary; better would be to substitute with values end -""" - get_limit_cycles(eom::HarmonicEquation, swept, fixed, ω_lc; kwargs...) - -Variant of `get_steady_states` for a limit cycle problem characterised by a Hopf frequency (usually called ω_lc) - -Solutions with ω_lc = 0 are labelled unphysical since this contradicts the assumption of distinct harmonic variables corresponding to distinct harmonics. -""" -function get_limit_cycles( - eom::HarmonicEquation, swept, fixed, ω_lc; explicit_Jacobian=false, kwargs... -) +function limit_cycle_problem(eom::HarmonicEquation, swept, fixed, ω_lc, explicit_Jacobian) prob = _cycle_Problem(eom, ω_lc) prob.jacobian = _gaugefixed_Jacobian( eom, @@ -103,20 +94,53 @@ function get_limit_cycles( sym_order=_free_symbols(prob, swept), rules=fixed, ) + return prob +end - result = get_steady_states( - prob, swept, fixed; method=:warmup, threading=true, classify_default=true, kwargs... - ) +""" + get_limit_cycles( + eom::HarmonicEquation, method::HarmonicBalanceMethod, swept, fixed, ω_lc; kwargs...) - _classify_limit_cycles!(result, ω_lc) +Variant of `get_steady_states` for a limit cycle problem characterised by a Hopf frequency (usually called ω_lc) - return result +Solutions with ω_lc = 0 are labelled unphysical since this contradicts the assumption of distinct harmonic variables corresponding to distinct harmonics. +""" +function get_limit_cycles( + eom::HarmonicEquation, swept, fixed, ω_lc; explicit_Jacobian=false, kwargs... +) + prob = limit_cycle_problem(eom, swept, fixed, ω_lc, explicit_Jacobian) + return get_limit_cycles(prob, WarmUp(), swept, fixed, ω_lc; kwargs...) end - function get_limit_cycles( - eom::HarmonicEquation, swept, fixed; limit_cycle_harmonic, kwargs... + eom::HarmonicEquation, + method::HarmonicBalanceMethod, + swept, + fixed, + ω_lc; + explicit_Jacobian=false, + kwargs..., ) - return get_limit_cycles(eom, swept, fixed, limit_cycle_harmonic; kwargs...) + prob = limit_cycle_problem(eom, swept, fixed, ω_lc, explicit_Jacobian) + return get_limit_cycles(prob, method, swept, fixed, ω_lc; kwargs...) +end +function get_limit_cycles(eom::HarmonicEquation, pairs::Dict, ω_lc; kwargs...) + swept = filter(x -> length(x[2]) > 1, pairs) + fixed = filter(x -> length(x[2]) == 1, pairs) + return get_limit_cycles(eom, swept, fixed, ω_lc; kwargs...) +end +function get_limit_cycles( + eom::HarmonicEquation, method::HarmonicBalanceMethod, pairs::Dict, ω_lc; kwargs... +) + swept = filter(x -> length(x[2]) > 1, pairs) + fixed = filter(x -> length(x[2]) == 1, pairs) + return get_limit_cycles(eom, method, swept, fixed, ω_lc; kwargs...) +end +function get_limit_cycles( + prob::Problem, method::HarmonicBalanceMethod, swept, fixed, ω_lc; kwargs... +) + result = get_steady_states(prob, method, swept, fixed; kwargs...) + _classify_limit_cycles!(result, ω_lc) + return result end # if abs(ω_lc) < tol, set all classifications to false diff --git a/src/LinearResponse/LinearResponse.jl b/src/LinearResponse/LinearResponse.jl index c3390fbe..4d8a6eb5 100644 --- a/src/LinearResponse/LinearResponse.jl +++ b/src/LinearResponse/LinearResponse.jl @@ -13,8 +13,18 @@ using LinearAlgebra: norm, eigvals, eigen, eigvecs using OrderedCollections: OrderedDict using HarmonicBalance +using HarmonicBalance: + Result, + HarmonicVariable, + HarmonicEquation, + DifferentialEquation, + StateDict, + get_variables, + get_independent_variables + using HarmonicBalance: var_name, + d, rearrange_standard, _remove_brackets, expand_derivatives, @@ -30,6 +40,7 @@ using HarmonicBalance: fourier_transform, declare_variable, is_rearranged + using ..HC_wrapper include("types.jl") diff --git a/src/LinearResponse/Lorentzian_spectrum.jl b/src/LinearResponse/Lorentzian_spectrum.jl index 44c1ebce..2fab36c7 100644 --- a/src/LinearResponse/Lorentzian_spectrum.jl +++ b/src/LinearResponse/Lorentzian_spectrum.jl @@ -80,9 +80,6 @@ function JacobianSpectrum(res::Result; index::Int, branch::Int, force=false) solution_dict = get_single_solution(res; branch=branch, index=index) λs, vs = eigen(res.jacobian(solution_dict)) - # convert OrderedDict to Dict - see Symbolics issue #601 - solution_dict = Dict(get_single_solution(res; index=index, branch=branch)) - for (j, λ) in enumerate(λs) eigvec = vs[:, j] # the eigenvector @@ -90,7 +87,7 @@ function JacobianSpectrum(res::Result; index::Int, branch::Int, force=false) for pair in _get_uv_pairs(hvars) u, v = hvars[pair] eigvec_2d = eigvec[pair] # fetch the relevant part of the Jacobian eigenvector - ωnum = real(unwrap(Symbolics.substitute(u.ω, solution_dict))) + ωnum = real(unwrap(Symbolics.substitute(u.ω, Dict(solution_dict)))) # ^ the harmonic (numerical now) associated to this harmonic variable # eigvec_2d is associated to a natural variable -> this variable gets Lorentzian peaks diff --git a/src/LinearResponse/jacobians.jl b/src/LinearResponse/jacobians.jl index 15e5ae15..706d8158 100644 --- a/src/LinearResponse/jacobians.jl +++ b/src/LinearResponse/jacobians.jl @@ -21,7 +21,7 @@ end " Obtain a Jacobian from a `DifferentialEquation` by first converting it into a `HarmonicEquation`. " function get_Jacobian(diff_eom::DifferentialEquation) - @variables T + Symbolics.@variables T harmonic_eq = get_harmonic_equations( diff_eom; slow_time=T, fast_time=first(get_independent_variables(diff_eom)) ) diff --git a/src/LinearResponse/response.jl b/src/LinearResponse/response.jl index 312dcb68..376d5743 100644 --- a/src/LinearResponse/response.jl +++ b/src/LinearResponse/response.jl @@ -7,7 +7,7 @@ This routine cannot accept a `HarmonicEquation` since there, some time-derivativ """ function get_response_matrix(diff_eq::DifferentialEquation, freq::Num; order=2)::Matrix - @variables T, i + Symbolics.@variables T, i time = get_independent_variables(diff_eq)[1] eom = harmonic_ansatz(diff_eq, time) @@ -33,7 +33,7 @@ Any substitution rules not specified in `res` can be supplied in `rules`." function ResponseMatrix(res::Result; rules=Dict()) # get the symbolic response matrix - @variables Δ + Symbolics.@variables Δ M = get_response_matrix(res.problem.eom.natural_equation, Δ; order=2) M = substitute_all(M, merge(res.fixed_parameters, rules)) symbols = _free_symbols(res) diff --git a/src/Problem.jl b/src/Problem.jl new file mode 100644 index 00000000..3d47115e --- /dev/null +++ b/src/Problem.jl @@ -0,0 +1,70 @@ + +""" +$(TYPEDEF) + +Holds a set of algebraic equations describing the steady state of a system. + +# Fields +$(TYPEDFIELDS) + +# Constructors +```julia +Problem(eom::HarmonicEquation; Jacobian=true) # find and store the symbolic Jacobian +Problem(eom::HarmonicEquation; Jacobian="implicit") # ignore the Jacobian for now, compute implicitly later +Problem(eom::HarmonicEquation; Jacobian=J) # use J as the Jacobian (a function that takes a Dict) +Problem(eom::HarmonicEquation; Jacobian=false) # ignore the Jacobian +``` +""" +mutable struct Problem + "The harmonic variables to be solved for." + variables::Vector{Num} + "All symbols which are not the harmonic variables." + parameters::Vector{Num} + "The input object for HomotopyContinuation.jl solver methods." + system::HC.System + "The Jacobian matrix (possibly symbolic). + If `false`, the Jacobian is ignored (may be calculated implicitly after solving)." + jacobian + "The HarmonicEquation object used to generate this `Problem`." + eom::HarmonicEquation + + function Problem(variables, parameters, system, jacobian) + return new(variables, parameters, system, jacobian) + end #incomplete initialization for user-defined symbolic systems + function Problem(variables, parameters, system, jacobian, eom) + return new(variables, parameters, system, jacobian, eom) + end +end + +Symbolics.get_variables(p::Problem)::Vector{Num} = get_variables(p.eom) + +function Base.show(io::IO, p::Problem) + println(io, length(p.system.expressions), " algebraic equations for steady states") + println(io, "Variables: ", join(string.(p.variables), ", ")) + println(io, "Parameters: ", join(string.(p.parameters), ", ")) + return println(io, "Symbolic Jacobian: ", !(p.jacobian == false)) +end + +# assume this order of variables in all compiled function (transform_solutions, Jacobians) +function _free_symbols(p::Problem, varied) + return cat(p.variables, collect(keys(OrderedDict(varied))); dims=1) +end + +""" Compile the Jacobian from `prob`, inserting `fixed_parameters`. + Returns a function that takes a dictionary of variables and `swept_parameters` to give the Jacobian.""" +function _compile_Jacobian( + prob::Problem, swept_parameters::ParameterRange, fixed_parameters::ParameterList +) + if prob.jacobian isa Matrix + compiled_J = compile_matrix( + prob.jacobian, _free_symbols(prob, swept_parameters); rules=fixed_parameters + ) + elseif prob.jacobian == "implicit" + compiled_J = LinearResponse.get_implicit_Jacobian( + prob, swept_parameters, fixed_parameters + ) # leave implicit Jacobian as is + else + return prob.jacobian + end + return compiled_J +end diff --git a/src/Result.jl b/src/Result.jl new file mode 100644 index 00000000..f9705f97 --- /dev/null +++ b/src/Result.jl @@ -0,0 +1,60 @@ +""" +$(TYPEDEF) + +Stores the steady states of a HarmonicEquation. + +# Fields +$(TYPEDFIELDS) + +""" +mutable struct Result + "The variable values of steady-state solutions." + solutions::Array{Vector{SteadyState}} + "Values of all parameters for all solutions." + swept_parameters::ParameterRange + "The parameters fixed throughout the solutions." + fixed_parameters::ParameterList + "The `Problem` used to generate this." + problem::Problem + "Maps strings such as \"stable\", \"physical\" etc to arrays of values, classifying the solutions (see method `classify_solutions!`)." + classes::Dict{String,Array} + "The Jacobian with `fixed_parameters` already substituted. Accepts a dictionary specifying the solution. + If problem.jacobian is a symbolic matrix, this holds a compiled function. + If problem.jacobian was `false`, this holds a function that rearranges the equations to find J + only after numerical values are inserted (preferable in cases where the symbolic J would be very large)." + jacobian::Function + "Seed used for the solver" + seed::UInt32 + + function Result(sol, swept, fixed, problem, classes, J, seed) + return new(sol, swept, fixed, problem, classes, J, seed) + end + Result(sol, swept, fixed, problem, classes) = new(sol, swept, fixed, problem, classes) + Result(sol, swept, fixed, problem) = new(sol, swept, fixed, problem, Dict([])) +end + +Symbolics.get_variables(res::Result)::Vector{Num} = get_variables(res.problem) + +function Base.show(io::IO, r::Result) + println(io, "A steady state result for ", length(r.solutions), " parameter points") + println(io, "\nSolution branches: ", length(r.solutions[1])) + println( + io, " of which real: ", sum(push!(any.(classify_branch(r, "physical")), false)) + ) + println( + io, " of which stable: ", sum(push!(any.(classify_branch(r, "stable")), false)) + ) + return println(io, "\nClasses: ", join(keys(r.classes), ", ")) +end + +# assume this order of variables in all compiled function (transform_solutions, Jacobians) +function _free_symbols(res::Result) + return cat(res.problem.variables, collect(keys(res.swept_parameters)); dims=1) +end + +# overload to use [] for indexing +Base.getindex(r::Result, idx::Int...) = get_single_solution(r, idx) +Base.size(r::Result) = size(r.solutions) + +branch_count(r::Result) = length(r.solutions[1]) +get_branch(r::Result, idx) = getindex.(r.solutions, idx) diff --git a/src/methods.jl b/src/methods.jl new file mode 100644 index 00000000..7d29b8d1 --- /dev/null +++ b/src/methods.jl @@ -0,0 +1,224 @@ +""" + HarmonicBalanceMethod + +Abstract type for harmonic balance methods. +""" +abstract type HarmonicBalanceMethod end + +""" + TotalDegree + +The Total Degree homotopy method. Performs a homotopy ``H(x, t) = γ t G(x) + (1-t) F(x)`` +from the trivial polynomial system ``xᵢ^{dᵢ} +aᵢ`` with the maximal degree ``dᵢ`` determined +by the [Bezout bound](https://en.wikipedia.org/wiki/B%C3%A9zout%27s_theorem). The method +guarantees to find all solutions, however, it comes with a high computational cost. See +[HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/guides/totaldegree/) +for more information. + +# Fields +$(TYPEDFIELDS) +""" +Base.@kwdef struct TotalDegree <: HarmonicBalanceMethod + """Complex multiplying factor of the start system G(x) for the homotopy""" + gamma::Complex = cis(2π * rand(Random.MersenneTwister(0xd8e5d8df))) + """Boolean indicating if threading is enabled.""" + thread::Bool = Threads.nthreads() > 1 + """Options for the tracker.""" + tracker_options::HomotopyContinuation.TrackerOptions = HomotopyContinuation.TrackerOptions() + """Options for the endgame.""" + endgame_options::HomotopyContinuation.EndgameOptions = HomotopyContinuation.EndgameOptions() + """Compilation options.""" + compile::Union{Bool,Symbol} = HomotopyContinuation.COMPILE_DEFAULT[] + """Seed for random number generation.""" + seed::UInt32 = 0xd8e5d8df +end + +""" + Polyhedral + +The Polyhedral homotopy method. This method constructs a homotopy based on the polyhedral +structure of the polynomial system. It is more efficient than the Total Degree method for +sparse systems, meaning most of the coefficients are zero. It can be especially useful if +you don't need to find the zero solutions (`only_non_zero = true`), resulting in speed up. +See [HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/guides/polyhedral/) +for more information. + +# Fields +$(TYPEDFIELDS) +""" +Base.@kwdef struct Polyhedral <: HarmonicBalanceMethod + """Boolean indicating if only non-zero solutions are considered.""" + only_non_zero::Bool = false + """Boolean indicating if threading is enabled.""" + thread::Bool = Threads.nthreads() > 1 + """Options for the tracker.""" + tracker_options::HomotopyContinuation.TrackerOptions = HomotopyContinuation.TrackerOptions() + """Options for the endgame.""" + endgame_options::HomotopyContinuation.EndgameOptions = HomotopyContinuation.EndgameOptions() + """Compilation options.""" + compile::Union{Bool,Symbol} = HomotopyContinuation.COMPILE_DEFAULT[] + """Seed for random number generation.""" + seed::UInt32 = 0xd8e5d8df +end + +""" + WarmUp + +The Warm Up method. This method prepares a warmup system using the parameter at `index` +perturbed by `perturbation_size` and performs a homotopy using the warmup system to all other +systems in the parameter sweep. It is very efficient for systems with less bifurcation in the +parameter sweep. The Warm Up method does not guarantee to find all solutions. See +[HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/guides/many-systems/) +for more information. + +# Fields +$(TYPEDFIELDS) +""" +Base.@kwdef struct WarmUp <: HarmonicBalanceMethod + """Size of the perturbation.""" + perturbation_size::ComplexF64 = 1e-6 + 1e-6 * im + """Index for the endpoint.""" + index::Union{Int,EndpointRanges.Endpoint} = EndpointRanges.iend ÷ 2 + """Boolean indicating if threading is enabled.""" + thread::Bool = Threads.nthreads() > 1 + """Options for the tracker.""" + tracker_options::HomotopyContinuation.TrackerOptions = HomotopyContinuation.TrackerOptions() + """Options for the endgame.""" + endgame_options::HomotopyContinuation.EndgameOptions = HomotopyContinuation.EndgameOptions() + """Compilation options.""" + compile::Union{Bool,Symbol} = HomotopyContinuation.COMPILE_DEFAULT[] + """Seed for random number generation.""" + seed::UInt32 = 0xd8e5d8df +end + +""" + thread(method::HarmonicBalanceMethod) -> Bool + +Returns whether threading is enabled for the given method. The number of available threads +is controlled by the environment variable `JULIA_NUM_THREADS`. +""" +thread(method::HarmonicBalanceMethod) = method.thread + +""" + tracker(method::HarmonicBalanceMethod) -> HomotopyContinuation.TrackerOptions + +Returns the tracker options for the given method. See `HomotopyContinuation.TrackerOptions` +for the available options. +""" +tracker(method::HarmonicBalanceMethod) = method.tracker_options + +""" + endgame(method::HarmonicBalanceMethod) -> HomotopyContinuation.EndgameOptions + +Returns the endgame options for the given method. See `HomotopyContinuation.EndgameOptions` +for the available options. +""" +endgame(method::HarmonicBalanceMethod) = method.endgame_options + +""" + compile(method::HarmonicBalanceMethod) -> Union{Bool,Symbol} + +Returns the compile options for the given method. If `true` then a system is compiled to a +straight line program for evaluation. This induces a compilation overhead. If `false` then +the generated program is only interpreted. This is slower than the compiled version, +but does not introduce compilation overhead. +""" +compile(method::HarmonicBalanceMethod) = method.compile + +""" + seed(method::HarmonicBalanceMethod) -> UInt32 + +Returns the seed for random number generation for the given method. +""" +seed(method::HarmonicBalanceMethod) = method.seed + +""" + alg_default_options(method::HarmonicBalanceMethod) -> NamedTuple + +Returns a named tuple of default algorithm options for the given method. +""" +function alg_default_options(method::HarmonicBalanceMethod) + return ( + threading=thread(method), + tracker_options=tracker(method), + endgame_options=endgame(method), + compile=compile(method), + seed=seed(method), + ) +end + +""" + alg_specific_options(method::TotalDegree) -> NamedTuple + +Returns a named tuple of specific algorithm options for the Total Degree method. +""" +alg_specific_options(method::TotalDegree) = (gamma=method.gamma,) + +""" + alg_specific_options(method::Polyhedral) -> NamedTuple + +Returns a named tuple of specific algorithm options for the Polyhedral method. +""" +alg_specific_options(method::Polyhedral) = (only_non_zero=method.only_non_zero,) + +""" + alg_specific_options(method::WarmUp) -> NamedTuple + +Returns a named tuple of specific algorithm options for the Warm Up method. +""" +function alg_specific_options(method::WarmUp) + return (perturbation_size=method.perturbation_size, index=method.index) +end + +""" + method_symbol(m::Polyhedral) -> Symbol + +Returns the symbol for the Polyhedral method identified by HomotopyContinuation/jl. +""" +method_symbol(m::Polyhedral) = :polyhedral + +""" + method_symbol(m::TotalDegree) -> Symbol + +Returns the symbol for the Total Degree method identified by HomotopyContinuation.jl. +""" +method_symbol(m::TotalDegree) = :total_degree + +""" + Base.show(io::IO, m::WarmUp) + +Displays information about the Warm Up method. +""" +function Base.show(io::IO, m::WarmUp) + println(io, "Warm up method:") + println(io, "perturbation_size: ", m.perturbation_size) + println(io, "Threading: ", thread(m)) + println(io, "Compile: ", compile(m)) + return println(io, "Seed: ", seed(m)) +end + +""" + Base.show(io::IO, m::TotalDegree) + +Displays information about the Total Degree method. +""" +function Base.show(io::IO, m::TotalDegree) + println(io, "Total degree method:") + println(io, "Gamma: ", m.gamma) + println(io, "Threading: ", thread(m)) + println(io, "Compile: ", compile(m)) + return println(io, "Seed: ", seed(m)) +end + +""" + Base.show(io::IO, m::Polyhedral) + +Displays information about the Polyhedral method. +""" +function Base.show(io::IO, m::Polyhedral) + println(io, "Polyhedral method:") + println(io, "Zero solutions: ", !m.only_non_zero) + println(io, "Threading: ", thread(m)) + println(io, "Compile: ", compile(m)) + return println(io, "Seed: ", seed(m)) +end diff --git a/src/precompilation.jl b/src/precompilation.jl index dd269803..36c98a1e 100644 --- a/src/precompilation.jl +++ b/src/precompilation.jl @@ -12,8 +12,7 @@ dEOM = DifferentialEquation(natural_equation + forces, x) add_harmonic!(dEOM, x, ω) harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); -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) -res = get_steady_states(prob, varied, fixed; show_progress=false) +res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false) diff --git a/src/solve_homotopy.jl b/src/solve_homotopy.jl index 2bed7e03..f5e198ff 100644 --- a/src/solve_homotopy.jl +++ b/src/solve_homotopy.jl @@ -1,25 +1,25 @@ """ - get_steady_states(prob::Problem, + get_steady_states(problem::HarmonicEquation, + method::HarmonicBalanceMethod, swept_parameters::ParameterRange, fixed_parameters::ParameterList; - method=:warmup, - threading = Threads.nthreads() > 1, show_progress=true, - sorting="nearest") + sorting="nearest", + classify_default=true) -Solves `prob` over the ranges specified by `swept_parameters`, keeping `fixed_parameters` constant. -`swept_parameters` accepts pairs mapping symbolic variables to arrays or `LinRange`. +Solves `problem` with the `method` over the ranges specified by `swept_parameters`, +keeping `fixed_parameters` constant. +`swept_parameters` accepts pairs mapping symbolic variables to arrays or ranges. `fixed_parameters` accepts pairs mapping symbolic variables to numbers. -Keyword arguments -- `method`: If `:warmup` (default), a problem similar to `prob` but with random complex parameters is first solved to find all non-singular paths. The subsequent tracking to find results for all `swept_parameters` is then much faster than the initial solving. If `method=:total_degree`, each parameter point is solved separately by tracking the maximum number of paths (employs a total degree homotopy). -This takes far longer but can be more reliable. -- `threading`: If `true`, multithreaded support is activated. The number of available threads is set by the environment variable `JULIA_NUM_THREADS`. -- `sorting`: the method used by `sort_solutions` to get continuous solutions branches. The current options are `"hilbert"` (1D sorting along a Hilbert curve), `"nearest"` (nearest-neighbor sorting) and `"none"`. +### Keyword arguments - `show_progress`: Indicate whether a progress bar should be displayed. +- `sorting`: the method used by `sort_solutions` to get continuous solutions branches. The current options are `"hilbert"` (1D sorting along a Hilbert curve), `"nearest"` (nearest-neighbor sorting) and `"none"`. +- `classify_default`: If `true`, the solutions will be classified using the default classification method. -Example: solving a simple harmonic oscillator ``m \\ddot{x} + γ \\dot{x} + ω_0^2 x = F \\cos(ωt)`` -to obtain the response as a function of ``ω`` +### Example +solving a simple harmonic oscillator +``m \\ddot{x} + γ \\dot{x} + ω_0^2 x = F \\cos(ωt)`` to obtain the response as a function of ``ω`` ```julia-repl # having obtained a Problem object, let's find steady states julia> range = ParameterRange(ω => LinRange(0.8,1.2,100) ) # 100 parameter sets to solve @@ -36,10 +36,10 @@ A steady state result for 100 parameter points ``` -It is also possible to create multi-dimensional solutions plots. +It is also possible to perfrom 2-dimensional sweeps. ```julia-repl # The swept parameters take precedence over fixed -> use the same fixed -julia> range = ParameterRange(ω => LinRange(0.8,1.2,100), F => LinRange(0.1,1.0,10) ) # 100x10 parameter sets +julia> range = ParameterRange(ω => range(0.8,1.2,100), F => range(0.1,1.0,10) ) # The swept parameters take precedence over fixed -> the F in fixed is now ignored julia> get_steady_states(problem, range, fixed) @@ -56,98 +56,77 @@ A steady state result for 1000 parameter points """ function get_steady_states( prob::Problem, + method::HarmonicBalanceMethod, swept_parameters::ParameterRange, fixed_parameters::ParameterList; - method=:warmup, - threading=Threads.nthreads() > 1, show_progress=true, sorting="nearest", classify_default=true, - seed=nothing, - kwargs..., ) - - # set seed if provided - !isnothing(seed) && Random.seed!(seed) + Random.seed!(seed(method)) # make sure the variables are in our namespace to make them accessible later declare_variable.(string.(cat(prob.parameters, prob.variables; dims=1))) - # prepare an array of vectors, each representing one set of input parameters - # an n-dimensional sweep uses an n-dimensional array - unique_fixed = filter_duplicate_parameters(swept_parameters, fixed_parameters) - variable_names = var_name.([keys(fixed_parameters)..., keys(swept_parameters)...]) any([var_name(var) ∈ variable_names for var in get_variables(prob)]) && error("Cannot fix one of the variables!") - input_array = _prepare_input_params(prob, swept_parameters, unique_fixed) - # feed the array into HomotopyContinuation, get back an similar array of solutions - raw = _get_raw_solution( - prob, - input_array; - sweep=swept_parameters, - method=method, - threading=threading, - show_progress=show_progress, - seed=seed, - kwargs..., + unique_fixed, input_array = _prepare_input_params( + prob, swept_parameters, fixed_parameters ) - - # extract all the information we need from results - #rounded_solutions = unique_points.(HomotopyContinuation.solutions.(getindex.(raw, 1)); metric = EuclideanNorm(), atol=1E-14, rtol=1E-8) - rounded_solutions = HC.solutions.(getindex.(raw, 1)) - all(isempty.(rounded_solutions)) ? error("No solutions found!") : nothing - solutions = pad_solutions(rounded_solutions) + solutions = get_solutions(prob, method, input_array; show_progress=show_progress) compiled_J = _compile_Jacobian(prob, swept_parameters, unique_fixed) - # a "raw" solution struct result = Result( - solutions, swept_parameters, unique_fixed, prob, Dict(), compiled_J, seed + solutions, swept_parameters, unique_fixed, prob, Dict(), compiled_J, seed(method) ) - # sort into branches if sorting != "no_sorting" sort_solutions!(result; sorting=sorting, show_progress=show_progress) - else - nothing end classify_default ? _classify_default!(result) : nothing return result end -function get_steady_states(p::Problem, swept, fixed; kwargs...) - return get_steady_states(p, ParameterRange(swept), ParameterList(fixed); kwargs...) -end -function get_steady_states(eom::HarmonicEquation, swept, fixed; kwargs...) - return get_steady_states(Problem(eom), swept, fixed; kwargs...) -end -function get_steady_states(p, pairs::Dict; kwargs...) +function get_steady_states( + p::Problem, method::HarmonicBalanceMethod, swept, fixed; kwargs... +) return get_steady_states( - p, - filter(x -> length(x[2]) > 1, pairs), - filter(x -> length(x[2]) == 1, pairs); - kwargs..., + p, method, ParameterRange(swept), ParameterList(fixed); kwargs... ) end - -""" Compile the Jacobian from `prob`, inserting `fixed_parameters`. - Returns a function that takes a dictionary of variables and `swept_parameters` to give the Jacobian.""" -function _compile_Jacobian( - prob::Problem, swept_parameters::ParameterRange, fixed_parameters::ParameterList +function get_steady_states( + eom::HarmonicEquation, method::HarmonicBalanceMethod, swept, fixed; kwargs... ) - if prob.jacobian isa Matrix - compiled_J = compile_matrix( - prob.jacobian, _free_symbols(prob, swept_parameters); rules=fixed_parameters - ) - elseif prob.jacobian == "implicit" - compiled_J = LinearResponse.get_implicit_Jacobian( - prob, swept_parameters, fixed_parameters - ) # leave implicit Jacobian as is + return get_steady_states(Problem(eom), method, swept, fixed; kwargs...) +end +function get_steady_states(eom::HarmonicEquation, pairs::Dict; kwargs...) + swept = filter(x -> length(x[2]) > 1, pairs) + fixed = filter(x -> length(x[2]) == 1, pairs) + return get_steady_states(eom, swept, fixed; kwargs...) +end +function get_steady_states( + eom::HarmonicEquation, method::HarmonicBalanceMethod, pairs::Dict; kwargs... +) + swept = filter(x -> length(x[2]) > 1, pairs) + fixed = filter(x -> length(x[2]) == 1, pairs) + return get_steady_states(eom, method, swept, fixed; kwargs...) +end +function get_steady_states(eom::HarmonicEquation, swept, fixed; kwargs...) + return get_steady_states(Problem(eom), WarmUp(), swept, fixed; kwargs...) +end + +function get_solutions(prob, method, input_array; show_progress) + raw = _get_raw_solution(prob, method, input_array; show_progress=show_progress) + + solutions = HC.solutions.(getindex.(raw, 1)) + if all(isempty.(solutions)) + @warn "No solutions found!" + return solutions else - return prob.jacobian + pad_solutions(solutions) end - return compiled_J end """ @@ -176,7 +155,7 @@ function _prepare_input_params( prob::Problem, sweeps::ParameterRange, fixed_parameters::ParameterList ) # sweeping takes precedence over fixed_parameters - fixed_parameters = filter_duplicate_parameters(sweeps, fixed_parameters) + unique_fixed = filter_duplicate_parameters(sweeps, fixed_parameters) # fixed order of parameters all_keys = cat(collect(keys(sweeps)), collect(keys(fixed_parameters)); dims=1) @@ -200,7 +179,7 @@ function _prepare_input_params( # order each parameter vector to match the order in prob input_array = getindex.(input_array, [permutation]) # HC wants arrays, not tuples - return tuple_to_vector.(input_array) + return unique_fixed, tuple_to_vector.(input_array) end "Remove occurrences of `sweeps` elements from `fixed_parameters`." @@ -213,77 +192,70 @@ function filter_duplicate_parameters(sweeps, fixed_parameters) end "A random warmup solution is computed to use as `start_parameters` in the homotopy." -function _solve_warmup(problem::Problem, params_1D, sweep; threading, show_progress) +function _solve_warmup(problem::Problem, method::WarmUp, params; show_progress) # complex perturbation of the warmup parameters - complex_pert = [1e-5 * randn(ComplexF64) for p in problem.parameters] - real_pert = ones(length(params_1D[1])) - warmup_parameters = params_1D[end ÷ 2] .* (real_pert + complex_pert) + options = alg_specific_options(method) + l = length(params[1]) + + perturbation = ones(l) + options[:perturbation_size] * randn(ComplexF64, l) + warmup_parameters = params[options[:index]] .* perturbation warmup_solution = HC.solve( problem.system; start_system=:total_degree, target_parameters=warmup_parameters, - threading=threading, show_progress=show_progress, + alg_default_options(method)..., ) return warmup_parameters, warmup_solution end "Uses HomotopyContinuation to solve `problem` at specified `parameter_values`." function _get_raw_solution( - problem::Problem, - parameter_values; - sweep=ParameterRange(), - method=:warmup, - threading=false, - show_progress=true, - seed=nothing, - kwargs..., + problem::Problem, method::WarmUp, parameter_values; show_progress ) - # HomotopyContinuation accepts 1D arrays of parameter sets - params_1D = reshape(parameter_values, :, 1) + warmup_parameters, warmup_solution = _solve_warmup( + problem, method, parameter_values; show_progress=show_progress + ) + result_full = HC.solve( + problem.system, + HC.solutions(warmup_solution); + start_parameters=warmup_parameters, + target_parameters=parameter_values, + show_progress=show_progress, + alg_default_options(method)..., + ) - if method == :warmup && !isempty(sweep) - warmup_parameters, warmup_solution = _solve_warmup( - problem, params_1D, sweep; threading=threading, show_progress=show_progress - ) - result_full = HC.solve( - problem.system, - HC.solutions(warmup_solution); - start_parameters=warmup_parameters, - target_parameters=parameter_values, - threading=threading, - show_progress=show_progress, - seed=seed, - kwargs..., + return reshape(result_full, size(parameter_values)...) +end + +"Uses HomotopyContinuation to solve `problem` at specified `parameter_values`." +function _get_raw_solution( + problem::Problem, method::Union{TotalDegree,Polyhedral}, parameter_values; show_progress +) + result_full = Array{Vector{Any},1}(undef, length(parameter_values)) + if show_progress + bar = Progress( + length(parameter_values); + dt=1, + desc="Solving via $method homotopy ...", + barlen=50, ) - elseif method == :total_degree || method == :polyhedral - result_full = Array{Vector{Any},1}(undef, length(parameter_values)) - if show_progress - bar = Progress( - length(parameter_values); - dt=1, - desc="Solving via $method homotopy ...", - barlen=50, - ) - end - for i in eachindex(parameter_values) # do NOT thread this - p = parameter_values[i] - show_progress ? ProgressMeter.next!(bar) : nothing - result_full[i] = [ - HC.solve( - problem.system; - start_system=method, - target_parameters=p, - threading=threading, - show_progress=false, - seed=seed, - ), - p, - ] - end - else - error("Unknown method: ", string(method)) + end + for i in eachindex(parameter_values) # do NOT thread this + p = parameter_values[i] + show_progress ? ProgressMeter.next!(bar) : nothing + result_full[i] = [ + HC.solve( + problem.system; + start_system=method_symbol(method), + target_parameters=p, + show_progress=false, + alg_default_options(method)..., + alg_specific_options(method)..., + ), + p, + ] end # reshape back to the original shape diff --git a/src/types.jl b/src/types.jl index 3defa00f..645a22ca 100644 --- a/src/types.jl +++ b/src/types.jl @@ -4,282 +4,6 @@ const StateDict = OrderedDict{Num,ComplexF64}; const SteadyState = Vector{ComplexF64}; const ParameterVector = Vector{Float64}; -""" -$(TYPEDEF) - -Holds differential equation(s) of motion and a set of harmonics to expand each variable. -This is the primary input for `HarmonicBalance.jl` ; after inputting the equations, the harmonics - ansatz needs to be specified using `add_harmonic!`. - -# Fields -$(TYPEDFIELDS) - -## Example -```julia-repl -julia> @variables t, x(t), y(t), ω0, ω, F, k; - -# equivalent ways to enter the simple harmonic oscillator -julia> DifferentialEquation(d(x,t,2) + ω0^2 * x - F * cos(ω*t), x); -julia> DifferentialEquation(d(x,t,2) + ω0^2 * x ~ F * cos(ω*t), x); - -# two coupled oscillators, one of them driven -julia> DifferentialEquation([d(x,t,2) + ω0^2 * x - k*y, d(y,t,2) + ω0^2 * y - k*x] .~ [F * cos(ω*t), 0], [x,y]); -``` -""" -mutable struct DifferentialEquation - """Assigns to each variable an equation of motion.""" - equations::OrderedDict{Num,Equation} - """Assigns to each variable a set of harmonics.""" - harmonics::OrderedDict{Num,OrderedSet{Num}} - - function DifferentialEquation(eqs) - return new(eqs, OrderedDict(var => OrderedSet() for var in keys(eqs))) - end - - # uses the above constructor if no harmonics defined - function DifferentialEquation(eqs::Vector{Equation}, vars::Vector{Num}) - return DifferentialEquation(OrderedDict(zip(vars, eqs))) - end - - # if expressions are entered instead of equations, automatically set them = 0 - function DifferentialEquation(exprs::Vector{Num}, vars::Vector{Num}) - return DifferentialEquation(exprs .~ Int(0), vars) - end - - function DifferentialEquation(eq::Equation, var::Num) - typerhs = typeof(eq.rhs) - typelhs = typeof(eq.lhs) - if eq.rhs isa AbstractVector || eq.lhs isa AbstractVector - throw( - ArgumentError( - "The equation is of the form $(typerhs)~$(typelhs) is not supported. Commenly one forgot to broadcast the equation symbol `~`.", - ), - ) - end - return DifferentialEquation([eq], [var]) - end - function DifferentialEquation(eq::Equation, vars::Vector{Num}) - typerhs = typeof(eq.rhs) - typelhs = typeof(eq.lhs) - throw( - ArgumentError( - "The variables are of type $(typeof(vars)). Whereas you gave one equation of type $(typerhs)~$(typelhs). Commenly one forgot to broadcast the equation symbol `~`.", - ), - ) - end - DifferentialEquation(lhs::Num, var::Num) = DifferentialEquation([lhs ~ Int(0)], [var]) -end - -function Base.show(io::IO, diff_eq::DifferentialEquation) - println(io, "System of ", length(keys(diff_eq.equations)), " differential equations") - println(io, "Variables: ", join(keys(diff_eq.equations), ", ")) - print(io, "Harmonic ansatz: ") - for var in keys(diff_eq.harmonics) - print(io, string(var), " => ", join(string.(diff_eq.harmonics[var]), ", ")) - print(io, "; ") - end - println(io, "\n") - return [println(io, eq) for eq in values(diff_eq.equations)] -end - -""" -$(TYPEDEF) - -Holds a variable stored under `symbol` describing the harmonic `ω` of `natural_variable`. - -# Fields -$(TYPEDFIELDS) -""" -mutable struct HarmonicVariable - """Symbol of the variable in the HarmonicBalance namespace.""" - symbol::Num - """Human-readable labels of the variable, used for plotting.""" - name::String - """Type of the variable (u or v for quadratures, a for a constant, Hopf for Hopf etc.)""" - type::String - """The harmonic being described.""" - ω::Num - """The natural variable whose harmonic is being described.""" - natural_variable::Num -end - -function Base.show(io::IO, hv::HarmonicVariable) - return println( - io, - "Harmonic variable ", - string.(hv.symbol) * " for harmonic ", - string(hv.ω), - " of ", - string(hv.natural_variable), - ) -end - -""" -$(TYPEDEF) - -Holds a set of algebraic equations governing the harmonics of a `DifferentialEquation`. - -# Fields -$(TYPEDFIELDS) -""" -mutable struct HarmonicEquation - """A set of equations governing the harmonics.""" - equations::Vector{Equation} - """A set of variables describing the harmonics.""" - variables::Vector{HarmonicVariable} - """The parameters of the equation set.""" - parameters::Vector{Num} - "The natural equation (before the harmonic ansatz was used)." - natural_equation::DifferentialEquation - - # use a self-referential constructor with _parameters - function HarmonicEquation(equations, variables, nat_eq) - return (x = new(equations, variables, Vector{Num}([]), nat_eq); - x.parameters = _parameters(x); - x) - end - function HarmonicEquation(equations, variables, parameters, natural_equation) - return new(equations, variables, parameters, natural_equation) - end -end - -function Base.show(io::IO, eom::HarmonicEquation) - println(io, "A set of ", length(eom.equations), " harmonic equations") - println(io, "Variables: ", join(string.(get_variables(eom)), ", ")) - println(io, "Parameters: ", join(string.(eom.parameters), ", ")) - println(io, "\nHarmonic ansatz: ", _show_ansatz(eom)) - println(io, "\nHarmonic equations:") - return [println(io, "\n", eq) for eq in eom.equations] -end - -"""Gives the relation between `var` and the underlying natural variable.""" -function _show_ansatz(var::HarmonicVariable) - t = var.natural_variable.val.arguments - t = length(t) == 1 ? string(t[1]) : error("more than 1 independent variable") - ω = string(var.ω) - terms = Dict("u" => "*cos(" * ω * t * ")", "v" => "*sin(" * ω * t * ")", "a" => "") - return string(string(var.symbol) * terms[var.type]) -end - -"""Gives the full harmonic ansatz used to construct `eom`.""" -function _show_ansatz(eom::HarmonicEquation) - output = "" - vars = unique(getfield.(eom.variables, :natural_variable)) - for nat_var in vars - # the Hopf variable (limit cycle frequency) does not contribute a term - harm_vars = filter( - x -> isequal(nat_var, x.natural_variable) && x.type !== "Hopf", eom.variables - ) - ansatz = join([_show_ansatz(var) for var in harm_vars], " + ") - output *= "\n" * string(nat_var) * " = " * ansatz - end - return output -end - -""" -$(TYPEDEF) - -Holds a set of algebraic equations describing the steady state of a system. - -# Fields -$(TYPEDFIELDS) - -# Constructors -```julia -Problem(eom::HarmonicEquation; Jacobian=true) # find and store the symbolic Jacobian -Problem(eom::HarmonicEquation; Jacobian="implicit") # ignore the Jacobian for now, compute implicitly later -Problem(eom::HarmonicEquation; Jacobian=J) # use J as the Jacobian (a function that takes a Dict) -Problem(eom::HarmonicEquation; Jacobian=false) # ignore the Jacobian -``` -""" -mutable struct Problem - "The harmonic variables to be solved for." - variables::Vector{Num} - "All symbols which are not the harmonic variables." - parameters::Vector{Num} - "The input object for HomotopyContinuation.jl solver methods." - system::HC.System - "The Jacobian matrix (possibly symbolic). - If `false`, the Jacobian is ignored (may be calculated implicitly after solving)." - jacobian - "The HarmonicEquation object used to generate this `Problem`." - eom::HarmonicEquation - - function Problem(variables, parameters, system, jacobian) - return new(variables, parameters, system, jacobian) - end #incomplete initialization for user-defined symbolic systems - function Problem(variables, parameters, system, jacobian, eom) - return new(variables, parameters, system, jacobian, eom) - end -end - -function Base.show(io::IO, p::Problem) - println(io, length(p.system.expressions), " algebraic equations for steady states") - println(io, "Variables: ", join(string.(p.variables), ", ")) - println(io, "Parameters: ", join(string.(p.parameters), ", ")) - return println(io, "Symbolic Jacobian: ", !(p.jacobian == false)) -end - -# assume this order of variables in all compiled function (transform_solutions, Jacobians) -function _free_symbols(p::Problem, varied) - return cat(p.variables, collect(keys(OrderedDict(varied))); dims=1) -end - -""" -$(TYPEDEF) - -Stores the steady states of a HarmonicEquation. - -# Fields -$(TYPEDFIELDS) - -""" -mutable struct Result - "The variable values of steady-state solutions." - solutions::Array{Vector{SteadyState}} - "Values of all parameters for all solutions." - swept_parameters::ParameterRange - "The parameters fixed throughout the solutions." - fixed_parameters::ParameterList - "The `Problem` used to generate this." - problem::Problem - "Maps strings such as \"stable\", \"physical\" etc to arrays of values, classifying the solutions (see method `classify_solutions!`)." - classes::Dict{String,Array} - "The Jacobian with `fixed_parameters` already substituted. Accepts a dictionary specifying the solution. - If problem.jacobian is a symbolic matrix, this holds a compiled function. - If problem.jacobian was `false`, this holds a function that rearranges the equations to find J - only after numerical values are inserted (preferable in cases where the symbolic J would be very large)." - jacobian::Function - "Seed used for the solver" - seed::Union{Nothing,UInt32} - - function Result(sol, swept, fixed, problem, classes, J, seed) - return new(sol, swept, fixed, problem, classes, J, seed) - end - Result(sol, swept, fixed, problem, classes) = new(sol, swept, fixed, problem, classes) - Result(sol, swept, fixed, problem) = new(sol, swept, fixed, problem, Dict([])) -end - -function Base.show(io::IO, r::Result) - println(io, "A steady state result for ", length(r.solutions), " parameter points") - println(io, "\nSolution branches: ", length(r.solutions[1])) - println(io, " of which real: ", sum(any.(classify_branch(r, "physical")))) - println(io, " of which stable: ", sum(any.(classify_branch(r, "stable")))) - return println(io, "\nClasses: ", join(keys(r.classes), ", ")) -end - -# assume this order of variables in all compiled function (transform_solutions, Jacobians) -function _free_symbols(res::Result) - return cat(res.problem.variables, collect(keys(res.swept_parameters)); dims=1) -end - -# overload to use [] for indexing -Base.getindex(r::Result, idx::Int...) = get_single_solution(r, idx) -Base.size(r::Result) = size(r.solutions) - -branch_count(r::Result) = length(r.solutions[1]) -get_branch(r::Result, idx) = getindex.(r.solutions, idx) - """ Represents a sweep of one or more parameters of a `HarmonicEquation`. diff --git a/test/API.jl b/test/API.jl index 02696e1e..f380b36c 100644 --- a/test/API.jl +++ b/test/API.jl @@ -14,18 +14,25 @@ using HarmonicBalance @test_throws ArgumentError DifferentialEquation(rhs ~ [F * cos(ω * t), 0], [x, y]) end -@testset begin - @variables ω1, ω2, ωₘ, t, ω, F, γ, λ, x(t), y(t) - eqs = [d(x, t, 2) + (ω1^2 - λ * cos(ωₘ * t)) * x + γ * d(x, t)] +@testset "get_steady_states API" begin + @variables ω1, t, ω, F, γ, λ, x(t), y(t) + eqs = [d(x, t, 2) + (ω1^2 - λ * cos(2 * ω * t)) * x + γ * d(x, t)] diff_eq = DifferentialEquation(eqs, [x]) add_harmonic!(diff_eq, x, ω) # drive frequency, close to ω1 harmonic_eq = get_harmonic_equations(diff_eq) + prob = HarmonicBalance.Problem(harmonic_eq) + varied = ω => range(0.7, 1.3, 100) - @test_throws MethodError get_steady_states(harmonic_eq, varied, threading=true) - @test_throws ArgumentError get_steady_states(harmonic_eq, Dict(varied), threading=true) + @test_throws MethodError get_steady_states(harmonic_eq, varied) + @test_throws ArgumentError get_steady_states(harmonic_eq, Dict(varied)) + + fixed = Dict(ω1 => 1.0, γ => 0.005, λ => 0.1) + @test_throws MethodError get_steady_states(prob, Dict()) + @test_throws MethodError get_steady_states(prob, varied, fixed) + r = get_steady_states(prob, HarmonicBalance.WarmUp(), varied, fixed) end @testset "forgot variable" begin diff --git a/test/ModelingToolkitExt.jl b/test/ModelingToolkitExt.jl index 690bc7b3..8f85bc6e 100644 --- a/test/ModelingToolkitExt.jl +++ b/test/ModelingToolkitExt.jl @@ -17,7 +17,7 @@ end ) fixed = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01) - param = ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) + param = HarmonicBalance.ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) sys = ODESystem(diff_eq) for p in string.([α, ω, ω0, F, γ]) @@ -33,8 +33,9 @@ end d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x ) + sys = ODESystem(harmonic_eq) fixed = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01) - param = ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) + param = HarmonicBalance.ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) ODEProblem(diff_eq, [1.0, 0.0], (0, 100), param) end @@ -50,7 +51,7 @@ end harmonic_eq = get_harmonic_equations(diff_eq) fixed = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01) - param = ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) + param = HarmonicBalance.ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) @testset "ODESystem" begin sys = ODESystem(harmonic_eq) diff --git a/test/hysteresis_sweep.jl b/test/hysteresis_sweep.jl index b6a8bb58..e49117f4 100644 --- a/test/hysteresis_sweep.jl +++ b/test/hysteresis_sweep.jl @@ -11,7 +11,8 @@ harmonic_eq = get_harmonic_equations(diff_eq) # implement ansatz to get harmonic fixed = (α => 1, ω0 => 1.0, γ => 0.005, F => 0.005, η => 0.2) # fixed parameters varied = ω => range(0.95, 1.1, 10) # range of parameter values -result = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) +method = HarmonicBalance.WarmUp(; seed=SEED) +result = get_steady_states(harmonic_eq, method, varied, fixed; show_progress=false) followed_branches, _ = follow_branch(1, result) followed_branches, _ = follow_branch(1, result; y="√(u1^2+v1^2)") diff --git a/test/krylov.jl b/test/krylov.jl index 4f1cae94..69813368 100644 --- a/test/krylov.jl +++ b/test/krylov.jl @@ -27,5 +27,5 @@ end fixed = (ω0 => 1.0, γ => 0.005, α => 1.0, η => 0, F => 0.0, ψ => 0.0, θ => 0.0) varied = (ω => range(0.99, 1.01, 5), λ => range(1e-6, 0.05, 5)) -res1 = get_steady_states(harmonic_eq1, varied, fixed; show_progress=false, seed=SEED); -res2 = get_steady_states(harmonic_eq2, varied, fixed; show_progress=false, seed=SEED); +res1 = get_steady_states(harmonic_eq1, varied, fixed; show_progress=false); +res2 = get_steady_states(harmonic_eq2, varied, fixed; show_progress=false); diff --git a/test/limit_cycle.jl b/test/limit_cycle.jl index 8d498db7..60950b68 100644 --- a/test/limit_cycle.jl +++ b/test/limit_cycle.jl @@ -15,10 +15,8 @@ import HarmonicBalance.LinearResponse.plot_linear_response fixed = () varied = μ => range(2, 3, 2) - - result = get_limit_cycles( - harmonic_eq, varied, fixed, ω_lc; show_progress=false, seed=SEED - ) + method = HarmonicBalance.WarmUp(; seed=SEED) + result = get_limit_cycles(harmonic_eq, method, varied, fixed, ω_lc; show_progress=false) @test sum(any.(classify_branch(result, "stable"))) == 4 @test sum(any.(classify_branch(result, "unique_cycle"))) == 1 @@ -51,5 +49,5 @@ end # varied = (ω => range(0.992, 0.995, 2)) # # results -# result = get_limit_cycles(harmonic_eq, varied, fixed, ω_lc; threading=true) +# result = get_limit_cycles(harmonic_eq, varied, fixed, ω_lc) # end diff --git a/test/linear_response.jl b/test/linear_response.jl index a862cb91..605cff62 100644 --- a/test/linear_response.jl +++ b/test/linear_response.jl @@ -12,7 +12,7 @@ harmonic_eq = get_harmonic_equations(diff_eq) fixed = (α => 1, ω0 => 1.0, γ => 1e-2, F => 1e-6) varied = ω => range(0.9, 1.1, 10) -result = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) +result = get_steady_states(harmonic_eq, varied, fixed; show_progress=false) @testset "first order" begin plot_linear_response( @@ -51,9 +51,7 @@ end fixed = (γ => 0.008, ω0 => 1.0, α => 1.0, F => 0.0, ωₚ => 1.0, λ => 0.016) varied = (ω => range(0.995, 1.005, 20)) - result_ω = get_steady_states( - harmonic_eq, varied, fixed; show_progress=false, seed=SEED - ) + result_ω = get_steady_states(harmonic_eq, varied, fixed; show_progress=false) @test_throws ErrorException plot_eigenvalues(result_ω; branch=1) end end diff --git a/test/methods.jl b/test/methods.jl new file mode 100644 index 00000000..f23053bf --- /dev/null +++ b/test/methods.jl @@ -0,0 +1,43 @@ +using HarmonicBalance +using Test + +@testset "WarmUp" begin + @variables α ω ω0 λ γ t x(t) + diff_eq = DifferentialEquation( + d(d(x, t), t) + γ * d(x, t) + ω0^2 * (1 - λ * cos(2 * ω * t)) * x + α * x^3, x + ) + + add_harmonic!(diff_eq, x, ω) # + harmonic_eq = get_harmonic_equations(diff_eq) + + fixed = HarmonicBalance.ParameterList((α => 1.0, ω0 => 1.1, λ => 0.01, γ => 0.01)) + varied = HarmonicBalance.ParameterRange((ω => range(0.9, 1.1, 20),)) + prob = HarmonicBalance.Problem(harmonic_eq) + + unique_fixed, input_array = HarmonicBalance._prepare_input_params(prob, varied, fixed) + @test length.(input_array) == fill(5, 20) + + method = WarmUp() + warmup_parameters, warmup_solution = HarmonicBalance._solve_warmup( + prob, method, input_array; show_progress=false + ) + diff = abs.((input_array[method.index] .- warmup_parameters)) + @test diff ≈ zeros(5) atol = real(method.perturbation_size) * 5 +end + +@testset "Polyhedral" begin + @variables α ω ω0 λ γ t x(t) + diff_eq = DifferentialEquation( + d(d(x, t), t) + γ * d(x, t) + ω0^2 * (1 - λ * cos(2 * ω * t)) * x + α * x^3, x + ) + + add_harmonic!(diff_eq, x, ω) # + harmonic_eq = get_harmonic_equations(diff_eq) + + fixed = HarmonicBalance.ParameterList((α => 1.0, ω0 => 1.1, λ => 0.001, γ => 0.01)) + varied = HarmonicBalance.ParameterRange((ω => range(0.9, 1.1, 10),)) + result = get_steady_states( + harmonic_eq, Polyhedral(; only_non_zero=true), varied, fixed; show_progress=false + ) + @test sum(reduce(vcat, classify_branch(result, "stable"))) == 0 +end diff --git a/test/parametron.jl b/test/parametron.jl index 7c5483a0..6167326c 100644 --- a/test/parametron.jl +++ b/test/parametron.jl @@ -16,6 +16,8 @@ dEOM = DifferentialEquation(natural_equation + forces, x) add_harmonic!(dEOM, x, ω) harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); +method = HarmonicBalance.WarmUp(; seed=SEED) + @testset "undriven parametron" begin fixed = (Ω => 1.0, γ => 1e-2, λ => 5e-2, F => 0, α => 1.0, η => 0.3, θ => 0, ψ => 0) varied = ω => range(0.9, 1.1, 20) @@ -30,11 +32,11 @@ harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); @test length(prob.variables) == 2 @test length(prob.parameters) == 9 @test length(prob.system.parameters) == 9 - res = get_steady_states(prob, varied, fixed; show_progress=false, seed=SEED) + res = get_steady_states(prob, method, varied, fixed; show_progress=false) end @testset "steady states" begin - res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) + res = get_steady_states(harmonic_eq, method, varied, fixed; show_progress=false) @test length(res.solutions) == 20 @test length(res.solutions[1]) == 5 @test sum(any.(classify_branch(res, "physical"))) == 5 @@ -53,13 +55,13 @@ harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); @testset "implicit jacobian" begin p = HarmonicBalance.Problem(harmonic_eq; Jacobian="implicit") - res = get_steady_states(p, varied, fixed; show_progress=false, seed=SEED) + res = get_steady_states(p, method, varied, fixed; show_progress=false) end @testset "2d sweep" begin # try to run a 2D calculation fixed = (Ω => 1.0, γ => 1e-2, F => 0, α => 1.0, η => 0.1, θ => 0, ψ => 0) varied = (ω => range(0.9, 1.1, 20), λ => range(0.01, 0.05, 20)) - res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) + res = get_steady_states(harmonic_eq, method, varied, fixed; show_progress=false) end end diff --git a/test/plotting.jl b/test/plotting.jl index 86135686..5388c098 100644 --- a/test/plotting.jl +++ b/test/plotting.jl @@ -17,7 +17,7 @@ harmonic_eq = get_harmonic_equations(dEOM); fixed = (ω0 => 1.0, γ => 1e-2, λ => 5e-2, α => 1.0, η => 0.3) varied = ω => range(0.9, 1.1, 10) -res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) +res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false) # plot 1D result plot(res; x="ω", y="u1"); @@ -27,7 +27,7 @@ plot_phase_diagram(res); fixed = (ω0 => 1.0, γ => 1e-2, α => 1.0, η => 0.3) varied = (ω => range(0.9, 1.1, 5), λ => range(0.01, 0.05, 5)) -res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) +res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false) # plot 2D result plot_phase_diagram(res); diff --git a/test/runtests.jl b/test/runtests.jl index 1081be93..2a835ac9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -61,6 +61,7 @@ end @testset "Computing steady states" begin include("parametron.jl") include("krylov.jl") + include("methods.jl") end @testset "Processing solutions" begin diff --git a/test/time_evolution.jl b/test/time_evolution.jl index 1369e3bd..125f9a70 100644 --- a/test/time_evolution.jl +++ b/test/time_evolution.jl @@ -15,11 +15,8 @@ forces = F * cos(ω * t + θ) dEOM = DifferentialEquation(natural_equation + forces, x) add_harmonic!(dEOM, x, ω) harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); -p = HarmonicBalance.Problem(harmonic_eq); fixed = (Ω => 1.0, γ => 1e-2, λ => 5e-2, F => 1e-3, α => 1.0, η => 0.3, θ => 0, ψ => 0) -varied = ω => range(0.9, 1.1, 10) -res = get_steady_states(p, varied, fixed; show_progress=false, seed=SEED) tspan = (0.0, 10) sweep = AdiabaticSweep(ω => (0.9, 1.1), tspan) # linearly interpolate between two values at two times diff --git a/test/transform_solutions.jl b/test/transform_solutions.jl index 0411d61b..bf68be3b 100644 --- a/test/transform_solutions.jl +++ b/test/transform_solutions.jl @@ -19,7 +19,7 @@ harmonic_eq = get_harmonic_equations(dEOM); fixed = (Ω => 1.0, γ => 1e-2, F => 1e-3, α => 1.0) varied = ω => range(0.9, 1.1, 10) -res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED); +res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false); transform_solutions(res, "u1^2+v1^2") transform_solutions(res, "√(u1^2+v1^2)"; realify=true)