Harmonic variable overhaul (#26)

* HarmonicVariable overhauled to hold a single variable (not u,v)

* Jacobian linear response docs

* Symbolics 4.4.1 capped

Author: jkosata

Project.toml

@@ -1,7 +1,7 @@
 name = "HarmonicBalance"
 uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 authors = ["Jan Kosata <[email protected]>", "Javier del Pino <[email protected]>"]
-version = "0.4.4"
+version = "0.5.0"
 
 [deps]
 BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"
@@ -42,7 +42,7 @@ ProgressMeter = "1.7.2"
 PyCall = "1.93.0"
 PyPlot = "2.10.0"
 SymbolicUtils = "0.19.7"
-Symbolics = "4.4.1"
+Symbolics = "= 4.4.1"
 julia = "1.7.0"
 
 [extras]

docs/src/background/stability_response.md

@@ -2,9 +2,17 @@
 
 The core of the harmonic balance method is expressing the system's behaviour in terms of Fourier components or _harmonics_. For an $N$-coordinate system, we _choose_ a set of $M_i$ harmonics to describe each coordinate $x_i$ : 
 ```math
+\begin{equation} \label{eq:ansatz}
 x_i(t) = \sum_{j=1}^{M_i} u_{i,j}  (T)  \cos(\omega_{i,j} t)+ v_{i,j} (T) \sin(\omega_{i,j} t) \:,
+\end{equation}
 ```
-This means the system is now described using a discrete set of variables $u_{i,j}$ and $v_{i,j}$. Constructing a vector $\mathbf{u}(T) = (u_{1,1}(T), v_{1,1}(T), \ldots u_{N,M_N}(T), v_{N, M_N}(T))$, we may obtain the _harmonic equations_ (see [an example of this procedure](@ref Duffing_harmeq))
+This means the system is now described using a discrete set of variables $u_{i,j}$ and $v_{i,j}$. Constructing the vector 
+```math
+\begin{equation} \label{eq:harmvar}
+\mathbf{u}(T) = (u_{1,1}(T), v_{1,1}(T), \ldots u_{N,M_N}(T), v_{N, M_N}(T))\,,
+\end{equation}
+``` 
+we may obtain the _harmonic equations_ (see [an example of this procedure](@ref Duffing_harmeq))
 ```math
 \begin{equation} \label{eq:harmeq}
 \frac{d\mathbf{u}(T)}{dT}  = \bar{\mathbf{F}} (\mathbf{u})
@@ -26,17 +34,84 @@ Eq. \eqref{eq:Jaceq} is exactly solvable for $\delta \mathbf{u}(T)$ given an ini
 
 ```math
 \begin{equation} \label{eq:fluct_evo}
-    \delta \mathbf{u}(T) = \sum_{r}(\mathbf{v}_r\cdot \delta\mathbf{u}(T_0))\hspace{1mm}\mathbf{v}_r e^{\lambda_r T}.
+    \delta \mathbf{u}(T) = \sum_{r} c_r \hspace{1mm}\mathbf{v}_r e^{\lambda_r T}.
 \end{equation}
 ```
 
 The dynamical behaviour near the steady states is thus governed by $e^{ \lambda_r T}$: if $\mathrm{Re}(\lambda_r)<0$ for all $\lambda_r$, the state $\mathbf{u}_0$ is stable. Conversely, if $\mathrm{Re}(\lambda_r)>0$ for at least one $\lambda_r$, the state is unstable - perturbations such as noise or a small applied drive will force the system away from $\mathbf{u}_0$. 
 
 
-### Linear response (WIP)
+### Linear response
+
+The response of a stable steady state to an additional oscillatory force, caused by weak probes or noise, is often of interest. It can be calculated by solving for the perturbation $\delta \mathbf{u}(T)$ in the presence of an additional drive term.
+
+```math
+\begin{equation} \label{eq:Jacforced}
+\frac{d}{dT} \left[\delta \mathbf{u}(T)\right] =  J(\mathbf{u}_0) \delta \mathbf{u}(T) + \boldsymbol{\xi} \,e^{i \Omega T}\,,
+\end{equation}
+```
+
+Suppose we have found an eigenvector of $J(\mathbf{u}_0)$ such that $J(\mathbf{u}) \mathbf{v} = \lambda \mathbf{v}$. To solve Eq. \eqref{eq:Jacforced}, we insert $\delta \mathbf{u}(T) = A(\Omega)\, \mathbf{v} e^{i \Omega T}$. Projecting each side of Eq. \eqref{eq:Jacforced} onto $\mathbf{v}$ gives
+
+```math
+A(\Omega) \left( i \Omega  - \lambda \right)  = \boldsymbol{\xi} \cdot \mathbf{v} \quad \implies \quad A(\Omega) = \frac{\boldsymbol{\xi} \cdot \mathbf{v} }{-\text{Re}[\lambda] + i \left( \Omega - \text{Im}[\lambda] \right)}
+```
+
+We see that each eigenvalue $\lambda$ results in a linear response that is a Lorentzian centered at $\Omega = \text{Im}[\lambda]$. Effectively, the linear response matches that of a harmonic oscillator with resonance frequency $\text{Im}[\lambda]$ and damping $\text{Re}[\lambda]$.
+
+Knowing the response of the harmonic variables $\mathbf{u}(T)$, what is the corresponding behaviour of the "natural" variables $x_i(t)$ ? To find this out, we insert the perturbation back into the harmonic ansatz \eqref{eq:ansatz}. Since we require real variables, let us use $\delta \mathbf{u}(T) = A(\Omega) \left( \mathbf{v} \, e^{i \Omega T} +   \mathbf{v}^* \, e^{-i \Omega T} \right)$. Plugging this into
+```math
+\begin{equation}
+\delta x_{i}(t) = \sum_{j=1}^{M_i} \delta{u}_{i,j}(t) \cos(\omega_{i,j} \,t) + \delta v_{i,j} (t) \sin(\omega_{i,j} \,t) 
+\end{equation}
+```  
+and multiplying out the sines and cosines gives
+```math
+\begin{align} \label{eq:deltax}
+\delta x_i(t) = \sum_{j=1}^{M_i} \bigg\{ \left( \text{Re}[\delta u_{i,j}] - \text{Im}[\delta v_{i,j}] \right) \,\cos[(\omega_{i,j} - \Omega) t]  \\
++ \left( \text{Im}[\delta u_{i,j}] + \text{Re}[\delta v_{i,j}] \right) \,\sin[(\omega_{i,j} - \Omega) t] \nonumber \\
++ \left( \text{Re}[\delta u_{i,j}] + \text{Im}[\delta v_{i,j}] \right) \,\cos[(\omega_{i,j} + \Omega) t] \nonumber \\
++ \left( -\text{Im}[\delta u_{i,j}] + \text{Re}[\delta v_{i,j}] \right) \,\sin[(\omega_{i,j} + \Omega) t] \bigg\} \nonumber
+\end{align}
+```
+where $\delta u_{i,j}$ and $\delta v_{i,j}$ are the components of $\delta \mathbf{u}$ corresponding to the respective harmonics $\omega_{i,j}$. 
+
+We see that a motion of the harmonic variables at frequency $\Omega$ appears as motion of $\delta x_i(t)$ at frequencies $\omega_{i,j}\pm \Omega$. 
+
+To make sense of this, we normalize the vector $\delta \mathbf{u}$ and use normalised components $\delta \hat{u}_{i,j}$ and $\delta \hat{v}_{i,j}$. We also define the Lorentzian distribution
+```math
+\begin{equation}
+L(x)_{x_0, \gamma} = \frac{1}{\left( x - x_0 \right)^2 + \gamma^2 }
+\end{equation}
+```
+We see that all components of $\delta x_i(t)$ [Eq. \eqref{eq:deltax}] are proportional to $L(\Omega)_{\text{Im}[\lambda], \text{Re}[\lambda]}$. The first and last two summands are Lorentzians centered at $\pm \Omega$ which oscillate at $\omega_{i,j} \pm \Omega$, respectively. From this, we can extract the linear response function in Fourier space, $\chi (\tilde{\omega})$
+
+```math
+\begin{multline}
+| \chi [\delta x _i](\tilde{\omega}) |^2 = \sum_{j=1}^{M_i} \bigg\{  \left[ \left( \text{Re}[\delta \hat{u}_{i,j}] - \text{Im}[\delta \hat{v}_{i,j}] \right)^2 + \left( \text{Im}[\delta \hat{u}_{i,j}] + \text{Re}[\delta \hat{v}_{i,j}] \right)^2 \right] L(\omega_{i,j} - \tilde{\omega})_{\text{Im}[\lambda], \text{Re}[\lambda]} \\
++ \left[ \left( \text{Re}[\delta \hat{u}_{i,j}] + \text{Im}[\delta \hat{v}_{i,j}] \right)^2 + \left(  \text{Re}[\delta \hat{v}_{i,j}] - \text{Im}[\delta \hat{u}_{i,j}] \right)^2 \right] L(\tilde{\omega} - \omega_{i,j})_{\text{Im}[\lambda], \text{Re}[\lambda]} \bigg\}
+\end{multline}
+```
+Keeping in mind that $L(x)_{x_0, \gamma} = L(x + \Delta)_{x_0 + \Delta, \gamma}$ and the normalization $\delta \hat{u}_{i,j}^2 + \delta \hat{v}_{i,j}^2 = 1$, we can rewrite this as
+```math
+\begin{equation}
+|\chi [\delta x _i](\tilde{\omega})|^2 = \sum_{j=1}^{M_i} \left( 1 + \alpha_{i,j} \right) L(\tilde{\omega})_{\omega_{i,j} - \text{Im}[\lambda], \text{Re}[\lambda]}
++ \left( 1 - \alpha_{i,j} \right) L(\tilde{\omega})_{\omega_{i,j} + \text{Im}[\lambda], \text{Re}[\lambda]}
+\end{equation}
+```
+
+where 
+```math
+\alpha_{i,j} = 2\left( \text{Im}[\delta \hat{u}_{i,j}] \text{Re}[\delta \hat{v}_{i,j}] - \text{Re}[\delta \hat{u}_{i,j}] \text{Im}[\delta \hat{v}_{i,j}] \right)
+```
+The above solution applies to every eigenvalue $\lambda$ of the Jacobian. It is now clear that the linear response function $\chi [\delta x _i](\tilde{\omega})$ contains for each eigenvalue $\lambda_r$ and harmonic $\omega_{i,j}$ : 
+
+* A Lorentzian centered at $\omega_{i,j}-\text{Im}[\lambda_r]$ with amplitude $1 + \alpha_{i,j}^{(r)}$ 
+* A Lorentzian centered at $\omega_{i,j}+\text{Im}[\lambda_r]$ with amplitude $1 - \alpha_{i,j}^{(r)}$ 
+
+_Sidenote:_ As $J$ a real matrix, there is an eigenvalue $\lambda_r^*$ for each $\lambda_r$. The maximum number of peaks in the linear response is thus equal to the dimensionality of $\mathbf{u}(T)$.
 
-The response of a stable steady state to an additional oscillatory force, caused by weak probes or noise, is often of interest. It can be calculated by solving for the perturbation $\delta \mathbf{u}$ in the presence of an additional drive term.
+The linear response of the system in the state $\mathbf{u}_0$ is thus fully specified by the complex eigenvalues and eigenvectors of $J(\mathbf{u}_0)$. In HarmonicBalance.jl, the module [LinearResponse](@ref linresp_man) creates a set of plottable [`Lorentzian`](@ref HarmonicBalance.LinearResponse.Lorentzian) objects to represent this.
 
-The linear response of the system in the state $\mathbf{u}_0$ is thus encoded in the complex eigenvalues and eigenvectors of $J(\mathbf{u}_0)$. 
 
 [Check out this example](@ref linresp_ex) of the linear response module of HarmonicBalance.jl

src/HarmonicEquation.jl

@@ -11,31 +11,36 @@ show(eom::HarmonicEquation) = show_fields(eom)
     harmonic_ansatz(eom::DifferentialEquation, time::Num; coordinates="Cartesian")
 
 Expand each variable of `diff_eom` using the harmonics assigned to it with `time` as the time variable.
-For each harmonic of each variable, an instance of `HarmonicVariable` (describing a pair of variables (u,v)) is automatically created and named. 
+For each harmonic of each variable, instance(s) of `HarmonicVariable` are automatically created and named. 
 
-`coordinates` allows for using different coordinate systems (e.g. 'polars') - CURRENTLY INACTIVE
 """
-harmonic_ansatz(eom::DifferentialEquation, time::Num; coordinates="Cartesian") = harmonic_ansatz(eom, time, Dict([[var, coordinates] for var in get_variables(eom)]))
-
-
-function harmonic_ansatz(diff_eom::DifferentialEquation, time::Num, trafo_types::Dict{Num,String})
+function harmonic_ansatz(diff_eom::DifferentialEquation, time::Num)
     !is_harmonic(diff_eom, time) && error("The differential equation is not harmonic in ", time, " !")
-    old_eqs = [diff_eom.equations[var] for var in get_variables(diff_eom)]
-    eqs, vars = deepcopy(old_eqs), []
-    rules = Dict()
-
-    idx_counter = 1 # keep count to label new variables
-    for var in get_variables(diff_eom)
-        to_substitute = 0 # holds all the subtitution rules
-        for ω in diff_eom.harmonics[var] # 2 new variables are created for each ω of each variable
-            these_names = coordinate_names[trafo_types[var]] # stores two strings denoting the variables as either "a", "ϕ", "u", "v"
-            unique_symbols = [  these_names[1] * string(idx_counter) , these_names[2] * string(idx_counter) ]
-            repl_rule, new_vars = _rotate_variable(var, ω, time, trafo_types[var], new_symbols=unique_symbols) # create the new variables
-            to_substitute += repl_rule
-            append!(vars, [new_vars]) 
-            idx_counter += 1
+    eqs = collect(values(diff_eom.equations))
+    rules, vars = Dict(), []
+
+    # keep count to label new variables
+    uv_idx = 1 
+    a_idx = 1
+
+    for nvar in get_variables(diff_eom) # sum over natural variables
+        to_substitute = 0 # combine all the subtitution rules for var
+        for ω in diff_eom.harmonics[nvar]
+            if !isequal(ω, 0) # nonzero harmonic - create u,v
+                rule_u, hvar_u = _create_harmonic_variable(nvar, ω, time, "u", new_symbol="u"*string(uv_idx))
+                rule_v, hvar_v = _create_harmonic_variable(nvar, ω, time, "v", new_symbol="v"*string(uv_idx))
+                rule = rule_u + rule_v
+                uv_idx += 1
+                append!(vars, [hvar_u, hvar_v])
+            else # zero harmonic - create a
+                rule, hvar = _create_harmonic_variable(nvar, ω, time, "a", new_symbol="a"*string(a_idx))
+                a_idx += 1
+                append!(vars, [hvar])
+            end
+            to_substitute += rule
         end
-        rules[var] = to_substitute
+        
+        rules[nvar] = to_substitute # total sub rule for nvar
     end
     eqs = substitute_all(eqs, rules)
     HarmonicEquation(eqs, Vector{HarmonicVariable}(vars), diff_eom)
@@ -155,15 +160,16 @@ $(TYPEDSIGNATURES)
 Return the independent variables (typically time) of `eom`.
 """
 function get_independent_variables(eom::HarmonicEquation)::Vector{Num}
-    dynamic_vars = flatten(getfield.(eom.variables, Symbol("symbols")))
+    dynamic_vars = flatten(getfield.(eom.variables, Symbol("symbol")))
     return flatten(unique([SymbolicUtils.arguments(var.val) for var in dynamic_vars]))
 end
 
 
 """
 $(TYPEDSIGNATURES)
 Extract the Fourier components of `eom` corresponding to the harmonics specified in `eom.variables`.
-For each harmonic of each variable, 2 equations are generated (cos and sin Fourier coefficients).
+For each non-zero harmonic of each variable, 2 equations are generated (cos and sin Fourier coefficients).
+For each zero (constant) harmonic, 1 equation is generated
 `time` does not appear in the resulting equations anymore.
 
 Underlying assumption: all time-dependences are harmonic.
@@ -176,20 +182,20 @@ end
 
 
 function fourier_transform!(eom::HarmonicEquation, time::Num)
-    avg_eqs = Vector{Equation}(undef, 2*length(eom.variables))
-
-    natural_variables = flatten([[x,x] for x in getfield.(eom.variables, :natural_variable)])
-    equation_assignments = flatten([findall(x -> isequal(x, var), collect(keys(eom.natural_equation.equations))) for var in natural_variables])
-    
-    for i in 1:length(avg_eqs)
-        j = Int(ceil(i/2))
-        ω = eom.variables[j].ω
-        eq = eom.equations[equation_assignments[i]]
-
-        if isodd(i)
-            avg_eqs[i] = fourier_cos_term(eq, ω, time)
-        else
-            avg_eqs[i] = fourier_sin_term(eq, ω, time)
+    avg_eqs = Vector{Equation}(undef, length(eom.variables))
+
+    # loop over the HarmonicVariables, each generates one equation
+    for (i, hvar) in enumerate(eom.variables)
+        # find the equation belonging to this variable
+        eq_idx = findall(x -> isequal(x, hvar.natural_variable), collect(keys(eom.natural_equation.equations)))
+        eq = eom.equations[eq_idx...]
+        # "type" is usually "u" or "v" (harmonic) or ["a"] (zero-harmonic)
+        if hvar.type == "u"
+            avg_eqs[i] = fourier_cos_term(eq, hvar.ω, time)
+        elseif hvar.type == "v" 
+            avg_eqs[i] = fourier_sin_term(eq, hvar.ω, time)
+        elseif hvar.type == "a"
+            avg_eqs[i] = fourier_cos_term(eq, 0, time) # pick out the constants
         end
     end
 

src/HarmonicVariable.jl

@@ -1,42 +1,35 @@
 import Symbolics.get_variables; export get_variables
 
 # pretty-printing
-display(vars::HarmonicVariable) = display(vars.names)
-display(vars::Vector{HarmonicVariable}) = display.(getfield.(vars, Symbol("names")))
+display(var::HarmonicVariable) = display(var.name)
+display(var::Vector{HarmonicVariable}) = display.(getfield.(var, Symbol("name")))
 
 
-# Each keyword gives two names for the new variables
-coordinate_names = Dict([["Cartesian", ["u", "v"]], ["polar", ["a", "ϕ"]]])
-
-
-function _coordinate_transform(new_vars, ω, t)
-    Dict([["Cartesian",  new_vars[1] * cos(ω*t) + new_vars[2] * sin(ω*t)],
-["polar", new_vars[1] * cos(ω*t + new_vars[2])]])
+function _coordinate_transform(new_var, ω, t, type)
+    coords = Dict([
+        "u" =>  new_var * cos(ω*t), 
+        "v" => new_var * sin(ω*t),
+        "a" => new_var])
+    return coords[type]
 end
 
 
-"Generates two variables describing the harmonic ω of a natural variable."
-function _rotate_variable(nat_var::Num, ω::Num, time::Num, transform::String; new_symbols::Vector{String})
-
-    types = coordinate_names[transform] # is the transformation "Cartesian" or "polar"
-    new_vars = [declare_variable(s, time) for s in new_symbols] # this holds the internal symbols
-
-    var_names_strings = [VDP_name * "_{" * var_name(nat_var) * "," * Base.replace(string(ω),"*"=>"") * "}" for VDP_name in types]
-    var_names = Dict(zip([Num(s) for s in new_symbols], var_names_strings))
-    repl_rule = _coordinate_transform(new_vars, ω, time)[transform] # a dictionary mapping the coordinate type to an expression
-    VDP = HarmonicVariable(new_vars, var_names, types, ω, nat_var)
-
-    repl_rule, VDP
+function _create_harmonic_variable(nat_var::Num, ω::Num, t::Num, type::String; new_symbol::String)
+    new_var = declare_variable(new_symbol, t) # this holds the internal symbol
+    name = type * "_{" * var_name(nat_var) * "," * Base.replace(string(ω),"*"=>"") * "}"
+    rule = _coordinate_transform(new_var, ω, t, type) # contribution of this harmonic variable to the natural variable
+    hvar = HarmonicVariable(new_var, name, type, ω, nat_var)
+    rule, hvar
 end
 
 
 ###
 # Functions for variable substutions and manipulation of HarmonicVariable
 ###
 
-function substitute_all(vars::HarmonicVariable, rules)
-    sym, freq = vars.symbols, vars.ω
-    HarmonicVariable(substitute_all(sym, rules), vars.names, vars.types, substitute_all(freq, rules), vars.natural_variable)
+function substitute_all(var::HarmonicVariable, rules)
+    sym, freq = var.symbol, var.ω
+    HarmonicVariable(substitute_all(sym, rules), var.name, var.type, substitute_all(freq, rules), var.natural_variable)
 end
 
 
@@ -46,7 +39,7 @@ substitute_all(vars::Vector{HarmonicVariable}, rules) = [substitute_all(var, rul
 "Returns the symbols of a `HarmonicVariable`."
 get_variables(vars::Vector{Num}) = unique(flatten([Num.(get_variables(x)) for x in vars]))
 
-get_variables(vars::HarmonicVariable) = unique(get_variables(vars.symbols))
+get_variables(var::HarmonicVariable) = Num.(get_variables(var.symbol))
 
 
 

src/Symbolics_utils.jl

@@ -241,7 +241,8 @@ function _fourier_term(x, ω, t, f)
     term = x * f(ω * t)
     term = trig_reduce(term)
     indep = get_independent(term, t)
-    2*Num(simplify_complex(Symbolics.expand(indep)))
+    ft = Num(simplify_complex(Symbolics.expand(indep)))
+    !isequal(ω, 0) ? 2*ft : ft # extra factor in case ω = 0 !
 end
 
 

src/modules/LinearResponse.jl

@@ -14,6 +14,7 @@ module LinearResponse
     import Base: *, show; export *, show
 
     include("LinearResponse/types.jl")
+    include("LinearResponse/utils.jl")
     include("LinearResponse/perturbations.jl")
     include("LinearResponse/jacobian_spectrum.jl")
     include("LinearResponse/linear_response.jl")