plot jacobian response of a parameter ODE sweep

src/HarmonicBalance.jl

@@ -57,6 +57,7 @@ module HarmonicBalance
 
     include("modules/LinearResponse.jl")
     using .LinearResponse
+    export plot_linear_response, plot_rotframe_jacobian_response
 
     include("modules/TimeEvolution.jl")
     using .TimeEvolution

src/modules/LinearResponse/plotting.jl

@@ -1,4 +1,5 @@
 using Plots, Latexify, ProgressMeter
+using HarmonicBalance: _set_Plots_default
 export plot_linear_response
 
 
@@ -19,6 +20,21 @@ function get_jacobian_response(res::Result, nat_var::Num, Ω_range, branch::Int;
     end
     C
 end
+function get_jacobian_response(res::Result, nat_var::Num, Ω_range, followed_branches::Vector{Int}; show_progress=true)
+
+    spectra = [JacobianSpectrum(res, branch=branch, index = i) for (i, branch) in pairs(followed_branches)]
+    C = Array{Float64, 2}(undef,  length(Ω_range), length(spectra))
+
+    if show_progress
+        bar = Progress(length(CartesianIndices(C)), 1, "Diagonalizing the Jacobian for each solution ... ", 50)
+    end
+    # evaluate the Jacobians for the different values of noise frequency Ω
+    for ij in CartesianIndices(C)
+        C[ij] = abs(evaluate(spectra[ij[2]][nat_var], Ω_range[ij[1]]))
+        show_progress ? next!(bar) : nothing
+    end
+    C
+end
 
 
 function get_linear_response(res::Result, nat_var::Num, Ω_range, branch::Int; order, show_progress=true)
@@ -65,7 +81,7 @@ function get_rotframe_jacobian_response(res::Result, Ω_range, branch::Int; show
             end
         end
         show_progress ? next!(bar) : nothing
-        
+
     end
     C
 end
@@ -91,8 +107,23 @@ X = Vector{Float64}(collect(values(res.swept_parameters))[1][stable])
 C = order == 1 ? get_jacobian_response(res, nat_var, Ω_range, branch, show_progress=show_progress) : get_linear_response(res, nat_var, Ω_range, branch; order=order, show_progress=show_progress)
 C = logscale ? log.(C) : C
 
-heatmap(X, Ω_range,  C; color=:viridis,
-    xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), HarmonicBalance._set_Plots_default..., kwargs...)
+xlabel = latexify(string(first(keys(res.swept_parameters)))); ylabel = latexify("Ω");
+heatmap(X, Ω_range,  C; color=:viridis, xlabel=xlabel, ylabel=ylabel, _set_Plots_default..., kwargs...)
+end
+function plot_linear_response(res::Result, nat_var::Num, followed_branches::Vector{Int}; Ω_range, logscale=false, show_progress=true, switch_axis = false, kwargs...)
+    length(size(res.solutions)) != 1 && error("1D plots of not-1D datasets are usually a bad idea.")
+
+    X = Vector{Float64}(collect(first(values(res.swept_parameters))))
+
+    C = get_jacobian_response(res, nat_var, Ω_range, followed_branches)
+    C = logscale ? log.(C) : C
+
+    xlabel = latexify(string(first(keys(res.swept_parameters)))); ylabel = latexify("Ω");
+    if switch_axis
+        heatmap(Ω_range, X, C'; color=:viridis, xlabel=ylabel, ylabel=xlabel, _set_Plots_default..., kwargs...)
+    else
+        heatmap(X, Ω_range, C; color=:viridis, xlabel=xlabel, ylabel=ylabel, _set_Plots_default..., kwargs...)
+    end
 end
 
 """
@@ -114,10 +145,10 @@ function plot_rotframe_jacobian_response(res::Result; Ω_range, branch::Int, log
     !isempty(findall(x->x==0, Ω_range)) && @warn("Probing with Ω=0 may lead to unexpected results")
 
     X = Vector{Float64}(collect(values(res.swept_parameters))[1][stable])
-    
+
     C = get_rotframe_jacobian_response(res, Ω_range, branch, show_progress=show_progress, damping_mod=damping_mod)
     C = logscale ? log.(C) : C
-    
+
     heatmap(X, Ω_range,  C; color=:viridis,
-       xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), HarmonicBalance._set_Plots_default..., kwargs...)
-end
+       xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), _set_Plots_default..., kwargs...)
+end

src/modules/LinearResponse/response.jl

@@ -1,7 +1,4 @@
 export get_response
-export plot_response
-
-
 
 """
     get_response_matrix(diff_eq::DifferentialEquation, freq::Num; order=2)
@@ -95,8 +92,3 @@ end
 # rotating frame into the lab frame
 _plusamp(uv) = norm(uv)^2 - 2*(imag(uv[1])*real(uv[2]) - real(uv[1])*imag(uv[2]))
 _minusamp(uv) = norm(uv)^2 + 2*(imag(uv[1])*real(uv[2]) - real(uv[1])*imag(uv[2]))
-
-
-
-
-

test/hysteresis_sweep.jl

@@ -11,8 +11,10 @@ fixed = (α => 1, ω0 => 1.0, γ => 0.005, F => 0.005, η => 0.2)   # fixed para
 varied = ω => range(0.95, 1.1, 10)           # range of parameter values
 result = get_steady_states(harmonic_eq, varied, fixed, show_progress=false)
 
-followed_branch, _ = follow_branch(1, result)
+followed_branches, _ = follow_branch(1, result)
 
-@test first(followed_branch) ≠ last(followed_branch)
+@test first(followed_branches) ≠ last(followed_branches)
 
 plot_1D_solutions_branch(1, result, x="ω", y="√(u1^2+v1^2)", show=false);
+
+plot_linear_response(result, x, followed_branches, Ω_range=range(0.9,1.1,10), show=false);

test/linear_response.jl

@@ -1,6 +1,4 @@
 using HarmonicBalance
-import HarmonicBalance.LinearResponse.plot_linear_response
-import HarmonicBalance.LinearResponse.plot_rotframe_jacobian_response
 
 @variables α, ω, ω0, F, γ, t, x(t);
 
@@ -14,4 +12,4 @@ result = get_steady_states(harmonic_eq, varied, fixed, show_progress=false)
 
 plot_linear_response(result, x, branch=1, Ω_range=range(0.9,1.1,10), order=1, logscale=true)
 
-plot_rotframe_jacobian_response(result, Ω_range=range(0.01,1,10), branch=1, logscale=true)
+plot_rotframe_jacobian_response(result, Ω_range=range(0.01,1.1,10), branch=1, logscale=true)