nonlinearoscillations/harmonicbalance.jl

A Julia package for solving nonlinear differential equations using the harmonic balance method.

HarmonicBalance.jl is a Julia package for solving nonlinear differential equations using the method of harmonic balance.

Installation

To install HarmonicBalance.jl, you can use the github repo or the Julia package manager,

using Pkg
Pkg.add("HarmonicBalance")

Documentation

For a detailed description of the package and examples, see the documentation.

This repo contains a collection of example notebooks.

Simple example

Let's find the steady states of a driven Duffing oscillator with nonlinear damping, its equation of motion is:

using HarmonicBalance
@variables α, ω, ω0, F, t, η, x(t) # declare constant variables and a function x(t)
diff_eq = DifferentialEquation(d(x,t,2) + ω0^2*x + α*x^3 + η*d(x,t)*x^2 ~ F*cos(ω*t), x)
add_harmonic!(diff_eq, x, ω) # specify the ansatz x = u(T) cos(ωt) + v(T) sin(ωt)

# implement ansatz to get harmonic equations
harmonic_eq = get_harmonic_equations(diff_eq)

# solve for a range of ω
result = get_steady_states(harmonic_eq, (ω => LinRange(0.9, 1.2, 100),
                                          α => 1., ω0 => 1.0, F => 0.01, η=>0.1))

A steady state result for 100 parameter points

Solution branches:   3
   of which real:    3
   of which stable:  2

Classes: stable, physical, Hopf, binary_labels

plot(result, "sqrt(u1^2 + v1^2)")

Citation

If you use HarmonicBalance.jl in your project, we kindly ask you to cite this paper, namely:

HarmonicBalance.jl: A Julia suite for nonlinear dynamics using harmonic balance Jan Košata, Javier del Pino, Toni L. Heugel, Oded Zilberberg SciPost Phys. Codebases 6 (2022)

harmonicbalance.jl's Issues

stable argument in plot_2D_phase_diagram_interactive ignored in the 1D cuts when cut_type="transform"

function get_harmonic_equations does not take into account hyperbolic sine

Hi,

While solving a set of differential equations:

free_eq = [m1*(d(x1,t,2)) + k1x1 + (k2/alpha)sinh(alpha(x1-x2)), m2(d(x2,t,2)) - (k2/alpha)sinh(alpha(x1-x2))]

The function get_harmonic_equations basically ignores the hyperbolic sine as we can see below:

A set of 4 harmonic equations
Variables: u1(T), v1(T), u2(T), v2(T)
Parameters: k1, ω, m1, k2, alpha, F, m2

Harmonic ansatz:
x1(t) = u1cos(ωt) + v1sin(ωt)
x2(t) = u2cos(ωt) + v2sin(ωt)

Harmonic equations:

k1u1(T) + (2//1)m1ωDifferential(T)(v1(T)) - F - m1*(ω^2)*u1(T) ~ 0

k1v1(T) - m1(ω^2)v1(T) - (2//1)m1ωDifferential(T)(u1(T)) ~ 0

(2//1)m2ωDifferential(T)(v2(T)) - m2(ω^2)*u2(T) ~ 0

(-2//1)m2ωDifferential(T)(u2(T)) - m2(ω^2)*v2(T) ~ 0

I had to linearize sinh to get the correct answer haha.

Add a progress bar

Right now, the only progress bar comes from HomotopyContinuation.jl, especially for parallel computations this can get messy. Replace this by a custom progress bar which shows the solution classifications and branch matching.

The user should choose the font

The package should let the use choose the font he wants to use. Otherwise it can cause problem when using HarmonicBalance in practice. For example, I use a lot of unicodes in my variable naming. However, the forced computer modern font does not recognise the unicodes in ploting such that I get plotting results seen in the figure.

The solutions would be to not touch the font setting of Plots.jl in the function

function _set_Plots_default()
    font = "Computer Modern"
    default(linewidth=2, legend=:outerright)
    default(fontfamily=font, titlefont=font , tickfont=font)
end

Still plot stable steady states full and unstable dashed when `class=physical`

Plotting when the keyword class=default and class=physical gives different results although the same solution are shown.
Indeed running

using HarmonicBalance, Plots

@variables t x(t) # symbolic variables
@variables ω ω0 Γ λ β η
natural_equation = d(d(x, t), t) + Γ*d(x, t) + ω0^2*(1-λ*cos(2*ω*t))*x + β*x^3 + η*d(x, t)*x^2
diff_eq = DifferentialEquation(natural_equation, x)

add_harmonic!(diff_eq, x, ω);
harmonic_eq = get_harmonic_equations(diff_eq);

varied_stability = (ω => LinRange(0.98, 1.03, 100))
fixed_stability = Dict(ω0 => 1.0, Γ => 0.01, β => 1.0, η => 0.3, λ => 0.03)
result_stability = get_steady_states(harmonic_eq, varied_stability, fixed_stability, threading=true, random_warmup = false)

plot_stability = plot(result_stability, x="ω", y="sqrt(u1^2 + v1^2)")

yields

whereas,

plot_stability = plot(result_stability, x="ω", y="sqrt(u1^2 + v1^2)", class = "physical")

Errors in "Simple example"

Hi,

my name is Patrick and I am a student in Oded's course on parametric phenomena at the University of Constance. We got the possibility to try this Julia package and I found two minor errors in the "Simple example" on the code overview page when installing and testing it.

1. solutions = get_steady_states(problem, swept, fixed) gives: UndefVarError: "problem not defined"
   -> problem should be changed to harmonic_eq

2. plot_1D_solutions(solutions, x="ω", y="sqrt(u1^2 + v1^2)"]); returns: "syntax: unexpected "]" in argument list"
   -> the closing bracket "]" should be removed

Fixing these errors will help future users to have a smoother start with the package (especially if they are not used to Julia as I was before). Thanks!

Best,
Patrick

update documentation for 2d plotting

2d plotting requires a keyword branch

update doc

Making plot recipes instead of multiple dispatch functions

It is recommended by the Plots.jl package to make plot recipes instead of multiple dispatch functions. Maybe we should explore this route. This would also easily resolve issue #59.

https://github.com/JuliaPlots/Plots.jl/tree/master/RecipesBase
https://docs.juliaplots.org/latest/recipes/

Parallelize symbolic operations

A macro @maybethread is ready for this, but some parts of SymbolicUtils.jl are currently not threadsafe.

Typical bottlenecks:

• fourier_term when getting the harmonic equations
• get_Jacobian for the symbolic Jacobian

Being able to export the precompilation results to a file

Hello!

I just wanted to write here an idea that would be extremely useful for accelerating simulation times. Once we have already compiled our system of differential equations, and once we can start changing parameters freely without having to wait again that long time, it would be pretty neat having the possibility to save to an offline file this result. this way, if on a different occasion we want to continue researching that new system after we had already closed Julia, we would not have to wait again for the compilation of the exact system. Specially for big systems, this could save a lot of time each day :)

Thanks and congratulations for this great package.

Pablo EA1FID

symbolic Jacobian fails for fractions

When an Equation contains fractions (SymbolicUtils.Div terms), it gets incorrectly classified as nonlinear by Symbolics.jl . As a result, we cannot rearrange it to get the symbolic Jacobian.

Change DataStructures.jl dependency with OrderedCollections.jl

Discussed in #70

^{Originally posted by oameye October 25, 2022}
Hey hey, I noticed that we use DataStructures.jl, the parent of OrderedCollections.jl. However, I think, we only use the OrderedDict type defined in OrderedCollections.jl. Hence, we only need the OrderedCollections.jl dependency. Although, opting for OrderedCollections.jl would not reduce are loading time by much. I think it is just cleaner to only load what you use. This is also suggested in the SCIML style guide.

Error for harmonics with decimal prefactors

Hi,

Some frequency driven systems exhibits a "Period-Doubling" behavior, where a response exists at half the excitation frequency.
Would it be possible to search for half frequency reponses ?

I tried :

add_harmonic!(system_eqs, I, [0.5ω,ω])

But it ends up with the following error :

ERROR: MethodError: no method matching isless(::ComplexF64, ::Int64)
Closest candidates are:
  isless(::Union{StatsBase.PValue, StatsBase.TestStat}, ::Real) at C:\Users\brsinquin\.julia\packages\StatsBase\XgjIN\src\statmodels.jl:90
  isless(::ArfLike, ::Integer) at C:\Users\brsinquin\.julia\packages\Arblib\PdCpA\src\predicates.jl:131
  isless(::Static.StaticInteger{X}, ::Real) where X at C:\Users\brsinquin\.julia\packages\Static\Ldb7F\src\Static.jl:455
  ...

when executing : harmonic_eqs = get_harmonic_equations(system_eqs)

ODEProblem method for type DifferentialEquation

To look at the time evolution for the type HarmonicEquation you define a method for the function ODEProblem which transform the HarmonicEquation into the DifferentialEquations.ODEProblem type.

It would be nice to have the same for the DifferentialEquation type.

the plotting namespace is not automatically imported

Plotting functions need to be called as HarmonicBalance.plotting_function from them to work

Interactive plotting is no longer a feature

Tests for plotting

As the recent update showed it could be handy to have some testing function for the plotting capabilities of HarmonicBalance. However, not sure how to test plotting other than testing if the plot funtion run.

Jacobian spectrum density plots (pcolormesh) have strange behaviour for missing stable solutions

In systems where no stable solutions are found, the missing values form gaps in the 1D solution plots, which is correct. However, one would expect a corresponding missing section in the 2D density plots of the Jacobian noise spectrum. This does not take place, as pcolormesh seems to fill the gaps with a continuous colored section, which is artificial and can be misleading. A solution would be to the spectrum in the noise spectrum with NaNs or Infs if no stable solutions have been found for the corresponding parameter value.

Loading older files fails if versions change

If a new branch is made, old .jld2 files stop working even if no relevant changes are made.

Error null_J not a function in this workspace.

Solution: define a dummy function null_J at loading as well as saving.

Distinguishing numerical zeroes from complex numbers

When plotting this minimal example for a amplitude sweep

using HarmonicBalance, Plots
using Plots

@variables t x(t) 
@parameters ω ω0 F0 Γ λ θ β ψ η
natural_equation = d(d(x, t), t) + Γ*d(x, t) + ω0^2*(1-λ*cos(2*ω*t + ψ))*x + β*x^3 + η*d(x, t)*x^2
forces = F0*cos(ω*t + θ)
diff_eq = DifferentialEquation(natural_equation ~ forces, x)

add_harmonic!(diff_eq, x, ω);
harmonic_eq = get_harmonic_equations(diff_eq);

varied = ω => LinRange(0.75, 1.25, 50)
fixed = Dict(ω0 => 1.0, Γ => 0.01, F0 => 0.0, β => 1.0, λ => 0.03, η => 0.0, θ => 0.0, ψ => 0.0)
result = get_steady_states(harmonic_eq, varied, fixed, random_warmup = false)
plot1 = plot(result, "sqrt(u1^2+v1^2)")

I get a warning that some values with non-negligible complex parts have been projected on the real axis. Although the comples part are negliable (powers to the -51). It would be nice if only warning were issued if the complex part higher than some threshold. Something like 1e-3 maybe.

First-order Harmonic Balance support

Can I start with an input of coupled first order differential equations instead of second-order ones?

Plot 1D cuts through 2D solutions

Especially for Hopf calculations, 2D plots are sometimes desirable (1 parameter + drive frequency) where each 1D cut gives different physical information.
----> implement something to avoid running two separate 1D calculations

Requires.jl

Requires.jl exports a macro, @require, that allows you to specify that some code is conditional on having both packages available.

Could be cool to make plotting only available if the user loads Plots.jl. This will make package precompile time faster. Although I admit it is maybe a bit overkill.

favicon for the harmonicBalance documentation website

A favicon for the harmonicBalance documentation website is missing

Transform solutions doesn't work for time-dependent solutions

DivideError: integer division error

Is the integer division error shown in picture below caused by the parameter settings or code flaw?

Add the function method `plot!(plot, data::Result)`

Add a method to the HarmonicBalance Plot function where you can add data to plot by the notation plot!(plot, data::Result) just as in the normal Plots.jl syntax.

Test performance by using Unityper

The SymbolicUtils package uses @compactified from Unityper.jl for the dispatch which yields improved speeds. The downside is that it reduced the readability of the code. Something, we should decide to use or not.

parallelize sorting functionality

Sorting the solution into branches is currently single-core, which becomes the bottleneck in Hopf calculations.

add keyword to the `plot` function to only show specific branches.

Would be nice to have a feature to only plot specific branches of the solutions. For example, the following would only plot branch 1 and 3 from the "physical" amplitude plot.

plot_amplitude = plot(result_amplitude, x="ω", y="sqrt(u1^2 + v1^2)", branches = [1,3])

documentation: "hilbert" or "nearest" default sorting

In the documentation it is stated that "hilbert" is the default sorting algortihm in get-steady_states function.

However, in the code it is "nearest"

Is it a typo?

Cannot use constant "harmonics"

Nonlinear terms which are of even power generate constant terms. In such cases, one must include constant displacement besides the usual harmonics.
This currently does not work - a constant displacement generates one harmonic variable instead of two and must be treated as a special case.

Linear response spectrum (Jacobian) assumes stable solutions always exist

Only the parameter extrema are used in the extent argument, but missing values (solution) are not taken into account. A quick fix is filling the matrix for the linear response spectra with NaNs when a solution is not found (in this way they will be ignored by imshow, while matrix dimensions are consistent).

Allow keyword arguments in the interface with DifferentialEquations.jl

Kwargs in ODEProblem(f::ODEFunction,u0,tspan,p=NullParameters();kwargs...), allow setting integrator options

Allow arbitrary keyword arguments in plotting functions to be passed to PyPlot/Matplotlib calls

`ODECompositeSolutions` undefined

Jacobian eigenvalue plotter is no longer supported in 0.60

Jacobian eigenvalues can still be accessed from get_Jacobian

Use precompilation

Just-in-time compilation times are getting pretty crazy (minutes). Most of the time, similar functions are used - get_harmonic_equations and get_steady_states. Precompilation could make this much more convenient.

Starting resource: https://julialang.org/blog/2021/01/precompile_tutorial/

Debloat HarmonicBalance.jl

Problem

The current implementation of HarmonicBalance.jl package is becoming bloated and difficult to maintain regarding building out the future NonlinearOscillations ecosystem. It would be beneficial to split it into multiple smaller packages in order to improve maintainability and ease of use. In that way, the packages under the ecosystem becomes more more lightweight, focused on a specific set of functionality.

Proposed solution

Mimic the SCIML ecosystem by creating a hierarchical structure of packages. Create a main interface package at the top of the hierarchy that serves as a low-dependency common interface for libraries downstream. Split the current package into smaller, focused packages that are designed to work well together. Another, maybe more relevant, ecosystem who does a similar thing would be QuantumOptics.jl.

Potential drawbacks

Implementation plan

• Create a new main interface package and move the common interface code into it. This would be the package now called HarmonicBalanceBase.jl.
• Identify the main functionalities of the current package that can be separated into smaller packages.
  • Computing slow-flow equations.
  • Solving slow-flow equations with Homotopy Continuation.
  • Time evolve and analyze slow-flow equations and the nonlinear oscillatory ODE's.
• Create new packages for each of the functionalities identified in the previous step and move the corresponding code into them.
• Find good names for all extra packages.
• Test the new packages to ensure they work well together and with the main interface package.

Additional information

• Here one can find some tips and tricks for making common interfaces between packages
• Package Reexport is commonly used.

Do allocate `ProgressBar` when not needed

Yes, In my tests only next! prints the ProgressBar. However, it is indeed not needed to allocate the ProgressBar. I will check if I can do that in a nice way this week. Maybe, it could also be nice to write some tests on the progress bar.

Originally posted by @oameye in #65 (comment)

load function fails for custom defined systems, time-independent systems

The structureeom is not created when calling Problem with a time-independent ODE (which is introduced manually, not derived symbolically). As a result, the load function fails,

`UndefRefError: access to undefined reference

Stacktrace:
[1] getproperty
@ ./Base.jl:42 [inlined]
[2] _parse_symbol_names(x::Problem)
@ HarmonicBalance ~/.julia/packages/HarmonicBalance/1VZoW/src/saving.jl:74
[3] _parse_symbol_names(x::Result)
@ HarmonicBalance ~/.julia/packages/HarmonicBalance/1VZoW/src/saving.jl:79
[4] load(filename::String)
@ HarmonicBalance ~/.julia/packages/HarmonicBalance/1VZoW/src/saving.jl:44`

Linear response spectrum in the rotating frame

Functions to retrieve the linear fluctuation spectrum in the lab frame already exist. Rotating frame spectrum can be useful to compare to e.g. lock-in amplifiers measurement

Sorting in 2d parameter space

Credit to @jdelpino:
The labeling of the branches differ for one dimensional cut through 2two dimensional results and just one dimensional sweep. Is there a way to label them the same?

The discrepancy is caused by your ambiguity when you sort in 2D parameter space. You can go along x, along y... and now we did it so it sorts according to the nearest neighbour. In 2D, you might be closer to a solution along the diagonal. (in parameter space) than to a solution that is in front of you along x. I think we could make the user decide if they want to sort along x or along y, instead of nearest neighbors.

Also the number of stable/unstable solutions differ?

This seems to be caused by the reordering, therefore the branch assignment, not doing a good job, since the number of solutions (curves) is probably the same in both cases.

Spaghetti plot function

I think with the switch form PyPlot to Plots we lost the plot_1D_solutions_spaghetti function in the process. Would be nice if feature would be added again. Here I share already a function I wrote for my notebooks. I will try later enhance and merge it with HarmonicBalance.jl. Note that all functions above the plot_1D_solutions_spaghetti function are already in place.

function _realify(x)
    !is_real(x) && !isnan(x) ? (@warn "Values with non-negligible complex parts have been projected on the real axis!", x) : nothing
    real(x)
end

_realify(v::Vector) = [_realify.(getindex.(v, k)) for k in 1:length(v[1])]
_realify(a::Array) = _realify.(a)

_vectorise(s::Vector) = s
_vectorise(s) = [s]

function _apply_mask(solns::Array{Vector{ComplexF64}},  booleans)
    factors = replace.(booleans, 0 => NaN)
    map(.*, solns, factors)
end
function _get_mask(res, classes, not_classes=[])
    classes == "all" && return fill(trues(length(res.solutions[1])), size(res.solutions))
    bools = vcat([res.classes[c] for c in _vectorise(classes)], [map(.!, res.classes[c]) for c in _vectorise(not_classes)])
    map(.*, bools...)
end

function plot_1D_solutions_spaghetti(res::Result)
    Z = res.swept_parameters[ω]
    u1 = transform_solutions(res,"u1")
    v1 = transform_solutions(res,"v1")
    v1_stable = _apply_mask(v1, _get_mask(res, ["physical", "stable"], [])) |> _realify
    u1_stable = _apply_mask(u1, _get_mask(res, ["physical", "stable"], [])) |> _realify
    v1_unstable = _apply_mask(v1, _get_mask(res, ["physical"], ["stable"])) |> _realify
    u1_unstable = _apply_mask(u1, _get_mask(res, ["physical"], ["stable"])) |> _realify
    p = plot(xlabel = "u1", ylabel = "v1"; _set_Plots_default...)
    for k in findall(x -> !all(isnan.(x)), v1_stable[1:end])
        plot!(p, u1_stable[k], v1_stable[k], Z, label = latexify(k))
    end
    for k in findall(x -> !all(isnan.(x)), v1_unstable[1:end])
        plot!(p, u1_unstable[k], v1_unstable[k], Z, label = latexify(k))
    end
    return p
end

`plot_2D_phase_diagram` fails if there is one solution only

Self-explanatory - when trying to plot a 2D phase diagram where all points have the same value, get
PyError
ValueError('array of sample points is empty')

Need of SnoopPrecompile for julia 1.9

I validated the compatibility with Julia v1.9.0-rc2. However, I noticed that the load and compile time of our package is actually slower with using SnoopPrecompile (26 sec) than without (22 sec). This is due to 1.9 improved compile speeds (biggest feature of the julia 1.9). SnoopPrecompile probably causes overhead.

A keyword to control the loading bar

It would be nice to opt out on showing a load bar as it is not so compatible with the jupyter notebook and you get sometimes get crazy long output,

The loading bar is a nice feature to have in the REPL tho 😄

`slow_flow!` incorrect for degrees higher than 2

Keeping higher-order time-derivatives in the harmonic may be useful but currently slow_flow! with degree > 2 leaves behind incorrectly transformed terms such as Differential(t)(Differential(T)(u1(T))))

`classify_solutions` slow for custom functions

The runtime for classify_solutions! scales badly for 2d sweeps (e.g. 100x100 sweeps).

Here I list some impovements:

• Make the for loop over the solutions threaded with the already implemented maybethread.
• When applying multiple classifiers at once, we can implement the functionality of applying all classifiers in the same iteration/'for loop".

Instability in computing steady states for certain parameter ranges

When plotting the stability diagram of the parametron with the following code,

using HarmonicBalance, ModelingToolkit, Plots

@variables t x(t) # symbolic variables
@parameters ω ω0 F0 Γ λ θ β ψ η
natural_equation = d(d(x, t), t) + Γ*d(x, t) + ω0^2*(1-λ*cos(2*ω*t + ψ))*x + β*x^3 + η*d(x, t)*x^2
forces = F0*cos(ω*t + θ)
diff_eq = DifferentialEquation(natural_equation ~ forces, x)

add_harmonic!(diff_eq, x, ω);
harmonic_eq = get_harmonic_equations(diff_eq);

varied_stability = (ω => LinRange(0.95, 1.05, 100), λ => LinRange(1e-6, 0.06, 100))
fixed_stability = Dict(ω0 => 1.0, Γ => 0.01, F0 => 0.0, β => 1.0, η => 0.3, θ => 0.0, ψ => 0.0)
result_stability = get_steady_states(harmonic_eq, varied_stability, fixed_stability);
plot_stability = plot_phase_diagram(result_stability, class= "stable", title = "η=$(fixed_stability[η])")

we find instability in computing steady states for certain parameter ranges.

This is caused by the fact random_warmup is set to true (probably) by default in get_steady_states. It seems the algorithm is not converging to all solutions, as it is doing the parameter homotopy from a system of polynomial that was adequate for a given parameter, but is not working well when you change the parameters of your system. The way that we implemented a fix is to start tracking the solutions from a trivial polynomial system which we know has the maximum number of solutions possible. This you do setting random_warmup = false in get steadystates.

It doesn't seem to be a problem we can solve, it might be more of an issue for the HomotopyContinuation.jl people.

A fix for now is setting random_warmup = false in get_steady_states function.