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:30
   # parameters solved: 50
       # paths tracked: 250






Solving for 2500 parameters...   3%|▋                   |  ETA: 0:00:30
   # parameters solved: 73
       # paths tracked: 365






Solving for 2500 parameters...   4%|▊                   |  ETA: 0:00:29
   # parameters solved: 100
       # paths tracked: 500






Solving for 2500 parameters...   5%|█                   |  ETA: 0:00:29
   # parameters solved: 124
       # paths tracked: 620






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






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






Solving for 2500 parameters...   8%|█▌                  |  ETA: 0:00:30
   # parameters solved: 190
       # paths tracked: 950






Solving for 2500 parameters...   9%|█▊                  |  ETA: 0:00:30
   # parameters solved: 213
       # paths tracked: 1065






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






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






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






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






Solving for 2500 parameters...  13%|██▋                 |  ETA: 0:00:29
   # parameters solved: 324
       # paths tracked: 1620






Solving for 2500 parameters...  14%|██▊                 |  ETA: 0:00:28
   # parameters solved: 351
       # paths tracked: 1755






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






Solving for 2500 parameters...  16%|███▏                |  ETA: 0:00:28
   # parameters solved: 392
       # paths tracked: 1960






Solving for 2500 parameters...  17%|███▎                |  ETA: 0:00:28
   # parameters solved: 413
       # paths tracked: 2065






Solving for 2500 parameters...  18%|███▌                |  ETA: 0:00:27
   # parameters solved: 441
       # paths tracked: 2205






Solving for 2500 parameters...  19%|███▊                |  ETA: 0:00:27
   # parameters solved: 464
       # paths tracked: 2320






Solving for 2500 parameters...  20%|███▉                |  ETA: 0:00:27
   # parameters solved: 489
       # paths tracked: 2445






Solving for 2500 parameters...  20%|████▏               |  ETA: 0:00:26
   # parameters solved: 511
       # paths tracked: 2555






Solving for 2500 parameters...  21%|████▎               |  ETA: 0:00:26
   # parameters solved: 534
       # paths tracked: 2670






Solving for 2500 parameters...  22%|████▌               |  ETA: 0:00:26
   # parameters solved: 556
       # paths tracked: 2780






Solving for 2500 parameters...  23%|████▋               |  ETA: 0:00:26
   # parameters solved: 575
       # paths tracked: 2875






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






Solving for 2500 parameters...  25%|█████               |  ETA: 0:00:25
   # parameters solved: 622
       # paths tracked: 3110






Solving for 2500 parameters...  26%|█████▏              |  ETA: 0:00:25
   # parameters solved: 641
       # paths tracked: 3205






Solving for 2500 parameters...  27%|█████▍              |  ETA: 0:00:24
   # parameters solved: 667
       # paths tracked: 3335






Solving for 2500 parameters...  28%|█████▌              |  ETA: 0:00:24
   # parameters solved: 691
       # paths tracked: 3455






Solving for 2500 parameters...  28%|█████▊              |  ETA: 0:00:24
   # parameters solved: 712
       # paths tracked: 3560






Solving for 2500 parameters...  29%|█████▉              |  ETA: 0:00:24
   # parameters solved: 734
       # paths tracked: 3670






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






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






Solving for 2500 parameters...  32%|██████▌             |  ETA: 0:00:23
   # parameters solved: 805
       # paths tracked: 4025






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






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






Solving for 2500 parameters...  35%|███████             |  ETA: 0:00:22
   # parameters solved: 873
       # paths tracked: 4365






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






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






Solving for 2500 parameters...  38%|███████▌            |  ETA: 0:00:21
   # parameters solved: 942
       # paths tracked: 4710






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






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






Solving for 2500 parameters...  41%|████████▏           |  ETA: 0:00:20
   # parameters solved: 1017
       # paths tracked: 5085






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






Solving for 2500 parameters...  43%|████████▌           |  ETA: 0:00:19
   # parameters solved: 1064
       # paths tracked: 5320






Solving for 2500 parameters...  44%|████████▊           |  ETA: 0:00:19
   # parameters solved: 1088
       # paths tracked: 5440






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






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






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






Solving for 2500 parameters...  47%|█████████▍          |  ETA: 0:00:18
   # parameters solved: 1174
       # paths tracked: 5870






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: 1215
       # paths tracked: 6075






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






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






Solving for 2500 parameters...  52%|██████████▎         |  ETA: 0:00:16
   # parameters solved: 1289
       # paths tracked: 6445






Solving for 2500 parameters...  53%|██████████▌         |  ETA: 0:00:16
   # parameters solved: 1315
       # paths tracked: 6575






Solving for 2500 parameters...  53%|██████████▋         |  ETA: 0:00:16
   # parameters solved: 1334
       # paths tracked: 6670






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






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






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






Solving for 2500 parameters...  57%|███████████▌        |  ETA: 0:00:14
   # parameters solved: 1431
       # paths tracked: 7155






Solving for 2500 parameters...  58%|███████████▋        |  ETA: 0:00:14
   # parameters solved: 1459
       # paths tracked: 7295






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






Solving for 2500 parameters...  60%|████████████        |  ETA: 0:00:13
   # parameters solved: 1505
       # paths tracked: 7525






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






Solving for 2500 parameters...  62%|████████████▌       |  ETA: 0:00:13
   # parameters solved: 1555
       # paths tracked: 7775






Solving for 2500 parameters...  63%|████████████▋       |  ETA: 0:00:12
   # parameters solved: 1574
       # paths tracked: 7870






Solving for 2500 parameters...  64%|████████████▊       |  ETA: 0:00:12
   # parameters solved: 1599
       # paths tracked: 7995






Solving for 2500 parameters...  65%|█████████████       |  ETA: 0:00:12
   # parameters solved: 1619
       # paths tracked: 8095






Solving for 2500 parameters...  66%|█████████████▏      |  ETA: 0:00:11
   # parameters solved: 1640
       # paths tracked: 8200






Solving for 2500 parameters...  67%|█████████████▎      |  ETA: 0:00:11
   # parameters solved: 1663
       # paths tracked: 8315






Solving for 2500 parameters...  68%|█████████████▌      |  ETA: 0:00:11
   # parameters solved: 1688
       # paths tracked: 8440






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






Solving for 2500 parameters...  70%|█████████████▉      |  ETA: 0:00:10
   # parameters solved: 1740
       # paths tracked: 8700






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






Solving for 2500 parameters...  71%|██████████████▎     |  ETA: 0:00:09
   # parameters solved: 1786
       # paths tracked: 8930






Solving for 2500 parameters...  72%|██████████████▌     |  ETA: 0:00:09
   # parameters solved: 1810
       # paths tracked: 9050






Solving for 2500 parameters...  73%|██████████████▋     |  ETA: 0:00:09
   # parameters solved: 1835
       # paths tracked: 9175






Solving for 2500 parameters...  74%|██████████████▉     |  ETA: 0:00:09
   # parameters solved: 1855
       # paths tracked: 9275






Solving for 2500 parameters...  75%|███████████████     |  ETA: 0:00:08
   # parameters solved: 1879
       # paths tracked: 9395






Solving for 2500 parameters...  76%|███████████████▎    |  ETA: 0:00:08
   # parameters solved: 1901
       # paths tracked: 9505






Solving for 2500 parameters...  77%|███████████████▍    |  ETA: 0:00:08
   # parameters solved: 1922
       # paths tracked: 9610






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






Solving for 2500 parameters...  79%|███████████████▊    |  ETA: 0:00:07
   # parameters solved: 1965
       # paths tracked: 9825






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






Solving for 2500 parameters...  80%|████████████████▏   |  ETA: 0:00:07
   # parameters solved: 2010
       # paths tracked: 10050






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






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






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






Solving for 2500 parameters...  84%|████████████████▊   |  ETA: 0:00:05
   # parameters solved: 2101
       # paths tracked: 10505






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






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: 2168
       # paths tracked: 10840






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






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






Solving for 2500 parameters...  89%|█████████████████▉  |  ETA: 0:00:04
   # parameters solved: 2232
       # paths tracked: 11160






Solving for 2500 parameters...  90%|██████████████████  |  ETA: 0:00:03
   # parameters solved: 2253
       # paths tracked: 11265






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






Solving for 2500 parameters...  92%|██████████████████▍ |  ETA: 0:00:03
   # parameters solved: 2303
       # paths tracked: 11515






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






Solving for 2500 parameters...  94%|██████████████████▊ |  ETA: 0:00:02
   # parameters solved: 2349
       # paths tracked: 11745






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






Solving for 2500 parameters...  96%|███████████████████▏|  ETA: 0:00:01
   # parameters solved: 2397
       # paths tracked: 11985






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






Solving for 2500 parameters...  98%|███████████████████▌|  ETA: 0:00:01
   # parameters solved: 2441
       # paths tracked: 12205






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






Solving for 2500 parameters...  99%|███████████████████▉|  ETA: 0:00:00
   # parameters solved: 2486
       # paths tracked: 12430






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.