docs: rm saving manuel (#401)

Author: oameye

docs/pages.jl

@@ -34,7 +34,6 @@ pages = [
         "Plotting" => "manual/plotting.md",
         "Time Evolution" => "manual/time_dependent.md",
         "Linear Response" => "manual/linear_response.md",
-        "Saving and Loading" => "manual/saving.md",
         "SciML Extension" => "manual/SciMLExt.md",
         ]
     ],

docs/src/examples/ab_initio_noise.md

@@ -9,6 +9,7 @@ This example demonstrates how to compute the spectra obtained from probing the s
 ````@example ab_initio_noise
 using HarmonicBalance, Plots
 using ModelingToolkit, StaticArrays, StochasticDiffEq, DSP
+using ModelingToolkit: setp
 ````
 
 We first define a  gelper function to compute power spectral density of the simulated response

docs/src/examples/cumulants_KPO.md

@@ -0,0 +1,104 @@
+```@meta
+EditURL = "../../../examples/cumulants_KPO.jl"
+```
+
+# Cumulant approximations for the Kerr parametric oscillator
+
+Let us compare the higher cumulant approximations for the Kerr parametric oscillator (KPO). The KPO is a model for a cavity with a Kerr nonlinearity and a parametric drive after one applied the Rotating Wave approximation. For this we use the `QuantumCumulants` and the `QuantumCumulantsExt` module in `HarmonicBalance`.
+
+````@example cumulants_KPO
+using QuantumCumulants, HarmonicBalance, Plots
+````
+
+## first order cumulant
+
+Let us define the Hamiltonian of the KPO using QuantumCumulants.
+
+````@example cumulants_KPO
+h = FockSpace(:cavity)
+@qnumbers a::Destroy(h)
+@variables Δ::Real U::Real G::Real κ::Real
+param = [Δ, U, G, κ]
+
+H_RWA = (-Δ + U) * a' * a + U * (a'^2 * a^2) / 2 - G * (a' * a' + a * a) / 2
+ops = [a, a']
+````
+
+We can now compute the mean-field equations for the KPO using the `meanfield` function from `QuantumCumulants`. We need the complete the system equations of moition using the `complete` function.
+
+````@example cumulants_KPO
+eqs_RWA = meanfield(ops, H_RWA, [a]; rates=[κ], order=1)
+eqs_completed_RWA = complete(eqs_RWA)
+````
+
+To compute the steady states of the KPO, we use `HarmonicBalance`. We define the HomotopyContinuation problem using the `Problem`.
+
+````@example cumulants_KPO
+fixed = (U => 0.001, κ => 0.002)
+varied = (Δ => range(-0.03, 0.03, 100), G => range(1e-5, 0.02, 100))
+problem_c1 = HarmonicBalance.Problem(eqs_completed_RWA, param, varied, fixed)
+````
+
+This gives us the phase
+
+````@example cumulants_KPO
+result = get_steady_states(problem_c1, WarmUp())
+plot_phase_diagram(result; class="stable")
+````
+
+## second order cumulant
+
+The next order cumulant can be computed by setting the `order` keyword.
+
+````@example cumulants_KPO
+ops = [a]
+eqs_RWA = meanfield(ops, H_RWA, [a]; rates=[κ], order=2)
+eqs_c2 = complete(eqs_RWA)
+problem_c2 = HarmonicBalance.Problem(eqs_c2, param, varied, fixed)
+````
+
+Which gives us the phase diagram
+
+````@example cumulants_KPO
+result = get_steady_states(problem_c2, WarmUp())
+plot_phase_diagram(result; class="stable", clim=(0, 4))
+````
+
+However, the phase diagram seems to wrong. Indeed, plotting, the photon number `a⁺aᵣ` shows that and additional state with  the negative photon number is stable. However, the this is unphysical.
+
+````@example cumulants_KPO
+fixed = (U => 0.001, κ => 0.002, G => 0.01)
+varied = (Δ => range(-0.03, 0.03, 200))
+problem_c2 = HarmonicBalance.Problem(eqs_c2, param, varied, fixed)
+result = get_steady_states(problem_c2, TotalDegree())
+plot(result; y="a⁺aᵣ", class="stable")
+````
+
+Let us classify the solutions having negative photon numbers. That way we can filter out these solutions.
+
+````@example cumulants_KPO
+classify_solutions!(result, "a⁺aᵣ < 0", "neg photon number");
+plot(result; y="aᵣ", class="stable", not_class="neg photon number")
+````
+
+## third order cumulant
+
+Similarly, for third order
+
+````@example cumulants_KPO
+ops = [a]
+eqs_RWA = meanfield(ops, H_RWA, [a]; rates=[κ], order=3)
+eqs_c3 = complete(eqs_RWA)
+
+fixed = (U => 0.001, κ => 0.002, G => 0.01)
+varied = (Δ => range(-0.03, 0.03, 50))
+problem_c3 = HarmonicBalance.Problem(eqs_c3, param, varied, fixed)
+result = get_steady_states(problem_c3, TotalDegree())
+classify_solutions!(result, "a⁺aᵣ < 0", "neg photon number");
+plot(result; y="a⁺aᵣ", class="stable")
+````
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+

docs/src/examples/harmonic_oscillator_KB_vs_HB.md

@@ -0,0 +1,92 @@
+```@meta
+EditURL = "../../../examples/harmonic_oscillator_KB_vs_HB.jl"
+```
+
+Harmonic oscillator: comparison of KB and HB methods
+
+````@example harmonic_oscillator_KB_vs_HB
+using HarmonicBalance, Plots
+
+@variables ω ω0 F γ t x(t)
+diff_eq = DifferentialEquation(d(x, t, 2) + ω0^2 * x + γ * d(x, t) ~ F * cos(ω * t), x)
+add_harmonic!(diff_eq, x, ω)
+````
+
+````@example harmonic_oscillator_KB_vs_HB
+krylov_eq1 = get_krylov_equations(diff_eq; order=1)
+````
+
+````@example harmonic_oscillator_KB_vs_HB
+krylov_eq2 = get_krylov_equations(diff_eq; order=2)
+````
+
+````@example harmonic_oscillator_KB_vs_HB
+harmonic_eq = get_harmonic_equations(diff_eq)
+harmonic_eq = rearrange_standard(harmonic_eq)
+````
+
+````@example harmonic_oscillator_KB_vs_HB
+varied = (ω => range(0.1, 1.9, 200)) # range of parameter values
+fixed = (ω0 => 1.0, γ => 0.05, F => 0.1) # fixed parameters
+result_krylov1 = get_steady_states(krylov_eq1, varied, fixed)
+result_krylov2 = get_steady_states(krylov_eq2, varied, fixed)
+result_harmonic = get_steady_states(harmonic_eq, varied, fixed);
+nothing #hide
+````
+
+````@example harmonic_oscillator_KB_vs_HB
+plot(
+    plot(result_krylov1; y="u1^2+v1^2", label="Krylov 1"),
+    plot(result_krylov2; y="u1^2+v1^2", label="Krylov 2"),
+    plot(result_harmonic; y="u1^2+v1^2", label="Harmonic");
+    layout=(3, 1),
+)
+````
+
+````@example harmonic_oscillator_KB_vs_HB
+plot(
+    plot(result_krylov1; y="u1", label="Krylov 1", legend=:best),
+    plot(result_krylov2; y="u1", label="Krylov 2", legend=:best),
+    plot(result_harmonic; y="u1", label="Harmonic", legend=:best);
+    layout=(3, 1),
+)
+````
+
+````@example harmonic_oscillator_KB_vs_HB
+plot(
+    plot(result_krylov1; y="v1", label="Krylov 1", legend=:best),
+    plot(result_krylov2; y="v1", label="Krylov 2", legend=:best),
+    plot(result_harmonic; y="v1", label="Harmonic", legend=:best);
+    layout=(3, 1),
+)
+````
+
+````@example harmonic_oscillator_KB_vs_HB
+plot(
+    plot_eigenvalues(result_krylov1, 1; title="Krylov 1", ylims=(-4, 4)),
+    plot_eigenvalues(result_krylov2, 1; title="Krylov 2", ylims=(-4, 4)),
+    plot_eigenvalues(result_harmonic, 1; title="Harmonic", ylims=(-4, 4));
+    layout=(3, 1),
+)
+````
+
+````@example harmonic_oscillator_KB_vs_HB
+plot(
+    plot_linear_response(
+        result_krylov1, x, 1; Ω_range=range(0.1, 1.9, 200), title="Krylov 1"
+    ),
+    plot_linear_response(
+        result_krylov2, x, 1; Ω_range=range(0.1, 1.9, 200), title="Krylov 2"
+    ),
+    plot_linear_response(
+        result_harmonic, x, 1; Ω_range=range(0.1, 1.9, 200), title="Harmonic"
+    );
+    layout=(3, 1),
+    clims=(0, 250),
+)
+````
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+

docs/src/examples/steady_state_sweep.md

@@ -0,0 +1,97 @@
+```@meta
+EditURL = "../../../examples/steady_state_sweep.jl"
+```
+
+# Steady state sweeps
+
+````@example steady_state_sweep
+using HarmonicBalance, SteadyStateDiffEq, ModelingToolkit
+using BenchmarkTools, Plots, StaticArrays, OrdinaryDiffEq, LinearAlgebra
+using HarmonicBalance: OrderedDict
+
+@variables α ω ω0 F γ η t x(t);
+
+diff_eq = DifferentialEquation(
+    d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) + η * d(x, t) * x^2 ~ F * cos(ω * t), x
+)
+add_harmonic!(diff_eq, x, ω)
+harmonic_eq = get_harmonic_equations(diff_eq)
+
+fixed = (ω0 => 1.0, γ => 1e-2, F => 0.02, α => 1.0, η => 0.3)
+ω_span = (0.8, 1.3);
+ω_range = range(ω_span..., 201);
+varied = ω => ω_range
+
+result_HB = get_steady_states(harmonic_eq, varied, fixed)
+plot(result_HB, "sqrt(u1^2+v1^2)")
+````
+
+## Steady state sweep using `SteadyStateDiffEq.jl`
+
+````@example steady_state_sweep
+param = OrderedDict(merge(Dict(fixed), Dict(ω => 1.1)))
+x0 = [0, 0.0];
+prob_ss = SteadyStateProblem(
+    harmonic_eq, x0, param; in_place=false, u0_constructor=x -> SVector(x...)
+)
+
+varied = 6 => ω_range
+result_ss = steady_state_sweep(
+    prob_ss, DynamicSS(Rodas5()); varied, abstol=1e-5, reltol=1e-5
+)
+
+plot(result_HB, "sqrt(u1^2+v1^2)")
+plot!(ω_range, norm.(result_ss); c=:gray, ls=:dash)
+````
+
+## Adiabatic evolution
+
+````@example steady_state_sweep
+timespan = (0.0, 50_000)
+sweep = AdiabaticSweep(ω => (0.8, 1.3), timespan) # linearly interpolate between two values at two times
+ode_problem = ODEProblem(harmonic_eq, fixed; u0=[0.01; 0.0], timespan, sweep)
+time_soln = solve(ode_problem, Tsit5(); saveat=250)
+
+plot(result_HB, "sqrt(u1^2+v1^2)")
+plot(time_soln.t, norm.(time_soln.u))
+````
+
+## using follow_branch
+
+````@example steady_state_sweep
+followed_branch, Ys = follow_branch(1, result_HB; y="√(u1^2+v1^2)")
+Y_followed_gr =
+    real.([Ys[param_idx][branch] for (param_idx, branch) in enumerate(followed_branch)]);
+
+plot(result_HB, "sqrt(u1^2+v1^2)")
+plot!(ω_range, Y_followed_gr; c=:gray, ls=:dash)
+````
+
+## comparison
+
+````@example steady_state_sweep
+@btime result_ss = steady_state_sweep(
+    prob_ss, DynamicSS(Rodas5()); varied, abstol=1e-5, reltol=1e-5
+)
+
+@btime time_soln = solve(ode_problem, Tsit5(); saveat=250)
+
+@btime begin
+    followed_branch, Ys = follow_branch(1, result_HB; y="√(u1^2+v1^2)")
+    Y_followed_gr =
+        real.([Ys[param_idx][branch] for (param_idx, branch) in enumerate(followed_branch)])
+end
+````
+
+Plotting them together gives:
+
+````@example steady_state_sweep
+plot(ω_range, norm.(result_ss))
+plot!(ω_range, norm.(time_soln.u))
+plot!(ω_range, Y_followed_gr)
+````
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+