Skip to content

Parametrically driven resonator

One of the most famous effects displaced by nonlinear oscillators is parametric resonance, where the frequency of the linear resonator is modulated in time Phys. Rev. E 94, 022201 (2016). In the following we analyse this system, governed by the equations

x¨(t)+γx˙(t)+Ω2(1λcos(2ωt+ψ))x+αx3+ηx2x˙+Fd(t)=0

where for completeness we also considered an external drive term Fd(t)=Fcos(ωt+θ) and a nonlinear damping term ηx2x˙

To implement this system in Harmonic Balance, we first import the library

julia
using HarmonicBalance, Plots

Subsequently, we type define parameters in the problem and the oscillating amplitude function x(t) using the variables macro from Symbolics.jl

julia
@variables ω₀ γ λ F η α ω t x(t)

natural_equation =
    d(d(x, t), t) +
    γ * d(x, t) +
    (ω₀^2 - λ * cos(2 * ω * t)) * x +
    α * x^3 +
    η * d(x, t) * x^2
forces = F * cos* t)
diff_eq = DifferentialEquation(natural_equation + forces, x)
System of 1 differential equations
Variables:       x(t)
Harmonic ansatz: x(t) => ;   

Differential(t)(Differential(t)(x(t))) + F*cos(t*ω) + Differential(t)(x(t))*γ + x(t)*(-cos(2t*ω)*λ + ω₀^2) + (x(t)^3)*α + (x(t)^2)*Differential(t)(x(t))*η ~ 0

Note that an equation of the form

mx¨+mω02(1λcos(2ωt+ψ))x+γx˙+αx3+ηx2x˙=Fcosωt

can be brought to dimensionless form by rescaling the units as described in Phys. Rev. E 94, 022201 (2016).

We are interested in studying the response of the oscillator to parametric driving and forcing. In particular, we focus on the first parametric resonance of the system, i.e. operating around twice the bare frequency of the undriven oscillator ω while the frequency of the external drive is also ω. For this purpose, we consider a harmonic ansatz which contains a single frequency: x(t)ucos(ωt)+vsin(ωt). In HarmonicBalance, we can do this via add_harmonic command:

julia
add_harmonic!(diff_eq, x, ω);

and replacing this by the time independent (averaged) equations of motion. This can be simply done by writing

julia
harmonic_eq = get_harmonic_equations(diff_eq)
A set of 2 harmonic equations
Variables: u1(T), v1(T)
Parameters: ω, α, γ, ω₀, λ, η, F

Harmonic ansatz: 
x(t) = u1(T)*cos(ωt) + v1(T)*sin(ωt)

Harmonic equations:

F - (1//2)*u1(T)*λ + (2//1)*Differential(T)(v1(T))*ω + Differential(T)(u1(T))*γ - u1(T)*(ω^2) + u1(T)*(ω₀^2) + v1(T)*γ*ω + (3//4)*(u1(T)^3)*α + (3//4)*(u1(T)^2)*Differential(T)(u1(T))*η + (1//2)*u1(T)*Differential(T)(v1(T))*v1(T)*η + (3//4)*u1(T)*(v1(T)^2)*α + (1//4)*(v1(T)^2)*Differential(T)(u1(T))*η + (1//4)*(u1(T)^2)*v1(T)*η*ω + (1//4)*(v1(T)^3)*η*ω ~ 0

Differential(T)(v1(T))*γ + (1//2)*v1(T)*λ - (2//1)*Differential(T)(u1(T))*ω - u1(T)*γ*ω - v1(T)*(ω^2) + v1(T)*(ω₀^2) + (1//4)*(u1(T)^2)*Differential(T)(v1(T))*η + (3//4)*(u1(T)^2)*v1(T)*α + (1//2)*u1(T)*v1(T)*Differential(T)(u1(T))*η + (3//4)*Differential(T)(v1(T))*(v1(T)^2)*η + (3//4)*(v1(T)^3)*α - (1//4)*(u1(T)^3)*η*ω - (1//4)*u1(T)*(v1(T)^2)*η*ω ~ 0

The output of these equations are consistent with the result found in the literature. Now we are interested in the linear response spectrum, which we can obtain from the solutions to the averaged equations (rotating frame) as a function of the external drive, after fixing all other parameters in the system. A call to get_steady_states then retrieves all steadystates found along the sweep employing the homotopy continuation method, which occurs in a complex space (see the nice HomotopyContinuation.jl docs)

1D parameters

We start with a varied set containing one parameter, ω,

julia
fixed = (ω₀ => 1.0, γ => 1e-2, λ => 5e-2, F => 1e-3, α => 1.0, η => 0.3)
varied = ω => range(0.9, 1.1, 100)

result = get_steady_states(harmonic_eq, varied, fixed)
A steady state result for 100 parameter points

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

Classes: stable, physical, Hopf

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 ω in sweep. This offers speed-up, but requires to be tested in each scenario against 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

julia
HarmonicBalance.save("parametron_result.jld2", result);

During the execution of get_steady_states, different solution branches are classified by their proximity in complex space, with subsequent filtering of real (physically acceptable solutions). In addition, the stability properties of each steady state is assessed from the eigenvalues of the Jacobian matrix. All this information can be succinctly represented in a 1D plot via

julia
plot(result; x="ω", y="sqrt(u1^2 + v1^2)")

The user can also introduce custom classes based on parameter conditions via classify_solutions!. Plots can be overlaid and use keywords from Plots, MarkdownAST.LineBreak()

julia
classify_solutions!(result, "sqrt(u1^2 + v1^2) > 0.1", "large")
plot(result, "sqrt(u1^2 + v1^2)"; class=["physical", "large"], style=:dash)
plot!(result, "sqrt(u1^2 + v1^2)"; not_class="large")

Alternatively, we may visualise all underlying solutions, including complex ones,

julia
plot(result, "sqrt(u1^2 + v1^2)"; class="all")

2D parameters

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 δ=ωLω0 and the parametric drive strength λ.

To perform a 2D sweep over driving frequency ω and parametric drive strength λ, we keep fixed from before but include 2 variables in varied

julia
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);

Solving for 2500 parameters...   2%|▍                   |  ETA: 0:00:33
   # parameters solved: 44
       # paths tracked: 220






Solving for 2500 parameters...   3%|▋                   |  ETA: 0:00:31
   # parameters solved: 72
       # paths tracked: 360






Solving for 2500 parameters...   4%|▊                   |  ETA: 0:00:30
   # parameters solved: 97
       # paths tracked: 485






Solving for 2500 parameters...   5%|█                   |  ETA: 0:00:30
   # parameters solved: 121
       # paths tracked: 605






Solving for 2500 parameters...   6%|█▏                  |  ETA: 0:00:30
   # parameters solved: 143
       # paths tracked: 715






Solving for 2500 parameters...   7%|█▍                  |  ETA: 0:00:30
   # parameters solved: 167
       # paths tracked: 835






Solving for 2500 parameters...   8%|█▌                  |  ETA: 0:00:29
   # parameters solved: 193
       # paths tracked: 965






Solving for 2500 parameters...   9%|█▊                  |  ETA: 0:00:29
   # parameters solved: 214
       # paths tracked: 1070






Solving for 2500 parameters...   9%|█▉                  |  ETA: 0:00:29
   # parameters solved: 237
       # paths tracked: 1185






Solving for 2500 parameters...  10%|██▏                 |  ETA: 0:00:29
   # parameters solved: 262
       # paths tracked: 1310






Solving for 2500 parameters...  11%|██▎                 |  ETA: 0:00:29
   # parameters solved: 285
       # paths tracked: 1425






Solving for 2500 parameters...  12%|██▍                 |  ETA: 0:00:29
   # parameters solved: 304
       # paths tracked: 1520






Solving for 2500 parameters...  13%|██▋                 |  ETA: 0:00:28
   # parameters solved: 328
       # paths tracked: 1640






Solving for 2500 parameters...  14%|██▉                 |  ETA: 0:00:28
   # parameters solved: 354
       # paths tracked: 1770






Solving for 2500 parameters...  15%|███                 |  ETA: 0:00:28
   # parameters solved: 376
       # paths tracked: 1880






Solving for 2500 parameters...  16%|███▎                |  ETA: 0:00:27
   # parameters solved: 403
       # paths tracked: 2015






Solving for 2500 parameters...  17%|███▍                |  ETA: 0:00:27
   # parameters solved: 428
       # paths tracked: 2140






Solving for 2500 parameters...  18%|███▋                |  ETA: 0:00:26
   # parameters solved: 452
       # paths tracked: 2260






Solving for 2500 parameters...  19%|███▊                |  ETA: 0:00:26
   # parameters solved: 472
       # paths tracked: 2360






Solving for 2500 parameters...  20%|████                |  ETA: 0:00:26
   # parameters solved: 497
       # paths tracked: 2485






Solving for 2500 parameters...  21%|████▏               |  ETA: 0:00:26
   # parameters solved: 514
       # paths tracked: 2570






Solving for 2500 parameters...  22%|████▍               |  ETA: 0:00:26
   # parameters solved: 540
       # paths tracked: 2700






Solving for 2500 parameters...  23%|████▌               |  ETA: 0:00:25
   # parameters solved: 569
       # paths tracked: 2845






Solving for 2500 parameters...  24%|████▊               |  ETA: 0:00:25
   # parameters solved: 592
       # paths tracked: 2960






Solving for 2500 parameters...  24%|████▉               |  ETA: 0:00:25
   # parameters solved: 612
       # paths tracked: 3060






Solving for 2500 parameters...  25%|█████▏              |  ETA: 0:00:24
   # parameters solved: 633
       # paths tracked: 3165






Solving for 2500 parameters...  26%|█████▎              |  ETA: 0:00:24
   # parameters solved: 656
       # paths tracked: 3280






Solving for 2500 parameters...  27%|█████▌              |  ETA: 0:00:24
   # parameters solved: 680
       # paths tracked: 3400






Solving for 2500 parameters...  28%|█████▋              |  ETA: 0:00:24
   # parameters solved: 704
       # paths tracked: 3520






Solving for 2500 parameters...  29%|█████▊              |  ETA: 0:00:23
   # parameters solved: 726
       # paths tracked: 3630






Solving for 2500 parameters...  30%|██████              |  ETA: 0:00:23
   # parameters solved: 751
       # paths tracked: 3755






Solving for 2500 parameters...  31%|██████▎             |  ETA: 0:00:23
   # parameters solved: 775
       # paths tracked: 3875






Solving for 2500 parameters...  32%|██████▌             |  ETA: 0:00:22
   # parameters solved: 806
       # paths tracked: 4030






Solving for 2500 parameters...  33%|██████▋             |  ETA: 0:00:22
   # parameters solved: 827
       # paths tracked: 4135






Solving for 2500 parameters...  34%|██████▊             |  ETA: 0:00:22
   # parameters solved: 847
       # paths tracked: 4235






Solving for 2500 parameters...  35%|███████             |  ETA: 0:00:21
   # parameters solved: 872
       # paths tracked: 4360






Solving for 2500 parameters...  36%|███████▏            |  ETA: 0:00:21
   # parameters solved: 894
       # paths tracked: 4470






Solving for 2500 parameters...  37%|███████▍            |  ETA: 0:00:21
   # parameters solved: 919
       # paths tracked: 4595






Solving for 2500 parameters...  38%|███████▌            |  ETA: 0:00:20
   # parameters solved: 940
       # paths tracked: 4700






Solving for 2500 parameters...  39%|███████▊            |  ETA: 0:00:20
   # parameters solved: 963
       # paths tracked: 4815






Solving for 2500 parameters...  39%|███████▉            |  ETA: 0:00:20
   # parameters solved: 986
       # paths tracked: 4930






Solving for 2500 parameters...  40%|████████            |  ETA: 0:00:20
   # parameters solved: 1007
       # paths tracked: 5035






Solving for 2500 parameters...  41%|████████▎           |  ETA: 0:00:19
   # parameters solved: 1030
       # paths tracked: 5150






Solving for 2500 parameters...  42%|████████▍           |  ETA: 0:00:19
   # parameters solved: 1054
       # paths tracked: 5270






Solving for 2500 parameters...  43%|████████▋           |  ETA: 0:00:19
   # parameters solved: 1080
       # paths tracked: 5400






Solving for 2500 parameters...  44%|████████▉           |  ETA: 0:00:18
   # parameters solved: 1103
       # paths tracked: 5515






Solving for 2500 parameters...  45%|█████████           |  ETA: 0:00:18
   # parameters solved: 1129
       # paths tracked: 5645






Solving for 2500 parameters...  46%|█████████▎          |  ETA: 0:00:18
   # parameters solved: 1156
       # paths tracked: 5780






Solving for 2500 parameters...  47%|█████████▍          |  ETA: 0:00:17
   # parameters solved: 1178
       # paths tracked: 5890






Solving for 2500 parameters...  48%|█████████▋          |  ETA: 0:00:17
   # parameters solved: 1196
       # paths tracked: 5980






Solving for 2500 parameters...  49%|█████████▊          |  ETA: 0:00:17
   # parameters solved: 1218
       # paths tracked: 6090






Solving for 2500 parameters...  50%|█████████▉          |  ETA: 0:00:17
   # parameters solved: 1242
       # paths tracked: 6210






Solving for 2500 parameters...  51%|██████████▏         |  ETA: 0:00:16
   # parameters solved: 1264
       # paths tracked: 6320






Solving for 2500 parameters...  51%|██████████▎         |  ETA: 0:00:16
   # parameters solved: 1286
       # paths tracked: 6430






Solving for 2500 parameters...  52%|██████████▌         |  ETA: 0:00:16
   # parameters solved: 1310
       # paths tracked: 6550






Solving for 2500 parameters...  53%|██████████▋         |  ETA: 0:00:15
   # parameters solved: 1332
       # paths tracked: 6660






Solving for 2500 parameters...  54%|██████████▉         |  ETA: 0:00:15
   # parameters solved: 1355
       # paths tracked: 6775






Solving for 2500 parameters...  55%|███████████         |  ETA: 0:00:15
   # parameters solved: 1377
       # paths tracked: 6885






Solving for 2500 parameters...  56%|███████████▎        |  ETA: 0:00:15
   # parameters solved: 1399
       # paths tracked: 6995






Solving for 2500 parameters...  57%|███████████▍        |  ETA: 0:00:14
   # parameters solved: 1420
       # paths tracked: 7100






Solving for 2500 parameters...  58%|███████████▌        |  ETA: 0:00:14
   # parameters solved: 1440
       # paths tracked: 7200






Solving for 2500 parameters...  59%|███████████▊        |  ETA: 0:00:14
   # parameters solved: 1463
       # paths tracked: 7315






Solving for 2500 parameters...  59%|███████████▉        |  ETA: 0:00:13
   # parameters solved: 1487
       # paths tracked: 7435






Solving for 2500 parameters...  61%|████████████▏       |  ETA: 0:00:13
   # parameters solved: 1513
       # paths tracked: 7565






Solving for 2500 parameters...  61%|████████████▎       |  ETA: 0:00:13
   # parameters solved: 1535
       # paths tracked: 7675






Solving for 2500 parameters...  62%|████████████▌       |  ETA: 0:00:12
   # parameters solved: 1558
       # paths tracked: 7790






Solving for 2500 parameters...  63%|████████████▊       |  ETA: 0:00:12
   # parameters solved: 1586
       # paths tracked: 7930






Solving for 2500 parameters...  64%|████████████▉       |  ETA: 0:00:12
   # parameters solved: 1610
       # paths tracked: 8050






Solving for 2500 parameters...  65%|█████████████       |  ETA: 0:00:11
   # parameters solved: 1632
       # paths tracked: 8160






Solving for 2500 parameters...  66%|█████████████▎      |  ETA: 0:00:11
   # parameters solved: 1654
       # paths tracked: 8270






Solving for 2500 parameters...  67%|█████████████▍      |  ETA: 0:00:11
   # parameters solved: 1677
       # paths tracked: 8385






Solving for 2500 parameters...  68%|█████████████▋      |  ETA: 0:00:11
   # parameters solved: 1703
       # paths tracked: 8515






Solving for 2500 parameters...  69%|█████████████▊      |  ETA: 0:00:10
   # parameters solved: 1725
       # paths tracked: 8625






Solving for 2500 parameters...  70%|██████████████      |  ETA: 0:00:10
   # parameters solved: 1750
       # paths tracked: 8750






Solving for 2500 parameters...  71%|██████████████▏     |  ETA: 0:00:10
   # parameters solved: 1771
       # paths tracked: 8855






Solving for 2500 parameters...  72%|██████████████▍     |  ETA: 0:00:09
   # parameters solved: 1796
       # paths tracked: 8980






Solving for 2500 parameters...  73%|██████████████▌     |  ETA: 0:00:09
   # parameters solved: 1820
       # paths tracked: 9100






Solving for 2500 parameters...  74%|██████████████▊     |  ETA: 0:00:09
   # parameters solved: 1841
       # paths tracked: 9205






Solving for 2500 parameters...  75%|██████████████▉     |  ETA: 0:00:08
   # parameters solved: 1864
       # paths tracked: 9320






Solving for 2500 parameters...  76%|███████████████▏    |  ETA: 0:00:08
   # parameters solved: 1890
       # paths tracked: 9450






Solving for 2500 parameters...  77%|███████████████▎    |  ETA: 0:00:08
   # parameters solved: 1914
       # paths tracked: 9570






Solving for 2500 parameters...  78%|███████████████▌    |  ETA: 0:00:07
   # parameters solved: 1938
       # paths tracked: 9690






Solving for 2500 parameters...  78%|███████████████▊    |  ETA: 0:00:07
   # parameters solved: 1961
       # paths tracked: 9805






Solving for 2500 parameters...  79%|███████████████▉    |  ETA: 0:00:07
   # parameters solved: 1985
       # paths tracked: 9925






Solving for 2500 parameters...  80%|████████████████▏   |  ETA: 0:00:06
   # parameters solved: 2011
       # paths tracked: 10055






Solving for 2500 parameters...  81%|████████████████▎   |  ETA: 0:00:06
   # parameters solved: 2034
       # paths tracked: 10170






Solving for 2500 parameters...  82%|████████████████▌   |  ETA: 0:00:06
   # parameters solved: 2059
       # paths tracked: 10295






Solving for 2500 parameters...  83%|████████████████▋   |  ETA: 0:00:06
   # parameters solved: 2081
       # paths tracked: 10405






Solving for 2500 parameters...  84%|████████████████▉   |  ETA: 0:00:05
   # parameters solved: 2105
       # paths tracked: 10525






Solving for 2500 parameters...  85%|█████████████████   |  ETA: 0:00:05
   # parameters solved: 2126
       # paths tracked: 10630






Solving for 2500 parameters...  86%|█████████████████▏  |  ETA: 0:00:05
   # parameters solved: 2146
       # paths tracked: 10730






Solving for 2500 parameters...  87%|█████████████████▍  |  ETA: 0:00:04
   # parameters solved: 2169
       # paths tracked: 10845






Solving for 2500 parameters...  88%|█████████████████▌  |  ETA: 0:00:04
   # parameters solved: 2191
       # paths tracked: 10955






Solving for 2500 parameters...  89%|█████████████████▊  |  ETA: 0:00:04
   # parameters solved: 2216
       # paths tracked: 11080






Solving for 2500 parameters...  90%|█████████████████▉  |  ETA: 0:00:03
   # parameters solved: 2240
       # paths tracked: 11200






Solving for 2500 parameters...  91%|██████████████████▏ |  ETA: 0:00:03
   # parameters solved: 2263
       # paths tracked: 11315






Solving for 2500 parameters...  91%|██████████████████▎ |  ETA: 0:00:03
   # parameters solved: 2283
       # paths tracked: 11415






Solving for 2500 parameters...  92%|██████████████████▌ |  ETA: 0:00:03
   # parameters solved: 2308
       # paths tracked: 11540






Solving for 2500 parameters...  93%|██████████████████▋ |  ETA: 0:00:02
   # parameters solved: 2332
       # paths tracked: 11660






Solving for 2500 parameters...  94%|██████████████████▉ |  ETA: 0:00:02
   # parameters solved: 2359
       # paths tracked: 11795






Solving for 2500 parameters...  95%|███████████████████ |  ETA: 0:00:02
   # parameters solved: 2382
       # paths tracked: 11910






Solving for 2500 parameters...  96%|███████████████████▎|  ETA: 0:00:01
   # parameters solved: 2403
       # paths tracked: 12015






Solving for 2500 parameters...  97%|███████████████████▍|  ETA: 0:00:01
   # parameters solved: 2427
       # paths tracked: 12135






Solving for 2500 parameters...  98%|███████████████████▋|  ETA: 0:00:01
   # parameters solved: 2447
       # paths tracked: 12235






Solving for 2500 parameters...  99%|███████████████████▊|  ETA: 0:00:00
   # parameters solved: 2472
       # paths tracked: 12360






Solving for 2500 parameters...  99%|████████████████████|  ETA: 0:00:00
   # parameters solved: 2497
       # paths tracked: 12485






Solving for 2500 parameters... 100%|████████████████████| Time: 0:00:33
   # parameters solved: 2500
       # paths tracked: 12500

Now, we count the number of solutions for each point and represent the corresponding phase diagram in parameter space. This is done using plot_phase_diagram. Only counting stable solutions,

julia
plot_phase_diagram(result_2D; class="stable")

In addition to phase diagrams, we can plot functions of the solution. The syntax is identical to 1D plotting. Let us overlay 2 branches into a single plot,

julia
# overlay branches with different colors
plot(result_2D, "sqrt(u1^2 + v1^2)"; branch=1, class="stable", camera=(60, -40))
plot!(result_2D, "sqrt(u1^2 + v1^2)"; branch=2, class="stable", color=:red)

Note that solutions are ordered in parameter space according to their closest neighbors. Plots can again be limited to a given class (e.g stable solutions only) through the keyword argument class.


This page was generated using Literate.jl.