Lab – Animating the Duffing Oscillator

Please write your name in this box

We encourage you to work together, and seek whatever help you need in completing this lab. Please list all people who you worked with or received help from

We encourage you to use any online resources that you find useful. Please list any URL’s that you consulted

We encourage you to use AI resources if you find them useful. If you made use of AI, please list the platform and the prompts which you made use of

Setup

In the last lab you explored the behavior of a chaotic driven nonlinear oscillator. Here you are going to use GLMakie to create a live animation that shows the phase space position of the oscillator as a function of time. It will have a Graphical User Interface (GUI) which will have a slider to adjust the frequency and amplitude of the drive. It will also have the ability to click on the graphic to set the phase-space position of the oscillator.

The Duffing oscillator is defined by

(93)\[\begin{align} \frac{dx}{dt}&=v\\ \frac{dv}{dt}&=-2\gamma v-\alpha x-\beta x^3 +F \cos(\omega t). \end{align}\]

Note, you will do a lot of copying and pasting in this Lab. I recommend either having a text file open to be a repository, or have two copies of this notebook open. To open a second copy in Jupyterlab, click on File and select New View for Notebook.

A second useful trick is the Simple toggle on the bottom left of Jupyterlab. If you click on it, your workspace is simplified to a single notebook. Click on it again, and it expands back to the setup that you had.

## Activity 1

Load GLMakie

using GLMakie

Our first step will be to create an object duffing which has the following attributes:

# Basic Scene info
fig # the Makie Figure object
ax # the Makie Axis object
# Time evolution info
u # the phase space position of the oscillator
t # the current time
dt # the step size
step # our stepper -- we will use rk4step as the default
dxdt # the derivitive rule
F # the amplitude of the force
ω # the frequency of the drive
α # parameter 
β # parameter 
γ # parameter 
# UI features
delay # how long to wait between frames
ball # the observable which will be plotted
traj # stores past points for drawing a "tail"
ptraj # not used in this lab
sliders # UI elements
run # button to turn on and off animation
isrunning # is simulation running?

We will make this a mutable struct – which means that we can change the entries after instantiation. The bare-bones way to do this is:

mutable struct duffing
    # Basic Scene info
    fig # the Makie Figure object
    ax # the Makie Axis object
    # Time evolution info
    xv # the phase space position of the oscillator
    t # the current time
    dt # the step size
    step # our stepper -- we will use rk4step as the default
    dxdt # the derivitive rule
    p # parameters of the differential equation -- a struct
      # with field F,ω,α,β,γ
    # UI features
    delay # how long to wait between frames
    ball # the observable which will be plotted
    traj # stores past points for drawing a "tail"
    ptraj # not used
    sliders # named tuple of UI elements to control force and frequency
            # Can again use the same struct that we used for p,
            # with attributes F,ω,α,β,γ
    run # button to turn on and off animation
    isrunning # is simulation running?
end

# Declaring the types of these attributes is optional
mutable struct dp
    F ::Float64
    ω ::Float64
    α ::Float64
    β ::Float64
    γ ::Float64
end

dp(;F,ω,α,β,γ)=dp(F,ω,α,β,γ)

dp

We need default values for the step and the dxdt. Here are the functions, coppied from previous work we did. I made one small change – setting them up so they take the parameters as an argument. This will be useful to us, as we will be changing the parameters on the fly.

function rk4step(x,p,dxdt,t,deltat)
    k1=dxdt(x,p,t)
    k2=dxdt(x+k1*(deltat/2),p,t+deltat/2)
    k3=dxdt(x+k2*(deltat/2),p,t+deltat/2)
    k4=dxdt(x+k3*deltat,p,t+deltat)
    return x+(k1+2*k2+2*k3+k4)*(deltat/6)
end

function duf_dxvdt(xv,p,t)
    x,v = xv
    return[v,-2*p.γ*v-p.α*x-p.β*x^3+p.F*cos(p.ω*t)]
end

duf_dxvdt (generic function with 1 method)

We now want to make a constructor, which is a function with the same name as our struct, which will create the duffing object, doing all of the necessary setup. Lets start with a bare-bones constructor, that just makes the scene and the ball object. We will just leave the other attributes as empty lists. We will make better constructors which fill things in as the lab progresses.

function duffing(;
        xv=[0.5,0.],
        p=dp(α=-1.,β=1.,γ=0.1,F=0.,ω=1.),
        t=0.,
        dt=0.01,
        delay=0.001,
        step=rk4step,
        dxdt=duf_dxvdt)
    # Create Figure object
    # 
    # Fill in details
    # must define `fig` and `ax`
    #
    xlims!(ax,-2,2) # set x axis range
    ylims!(ax,-2,2) # set y axis range
    # Create ball object, and add it to scene
    #
    # Fill in details
    # must define `ball`
    # this should be a Makie `observable`
    # it can be added to plot with `scatter!`
    # set its location to be the same as `xv`
    #
    d= duffing(
        fig #= the Makie Figure object that you created=#,
        ax #= the Makie Axis object that you created=#,
        xv #= initial conditions passed from constructor=#,
        t #= the current time passed from constructor=#,
        dt #= the step size passed from constructor=#,
        step #= our stepper passed from constructor=#,
        dxdt #= the derivitive rule passed from constructor=#,
        p #= passed from constructor=#,
        delay #= passed from constructor=#,
        ball #= the observable that you created =#,
        [] #=traj -- empty list =#,
        [] #=ptraj -- empty list =#,
        [] #=sliders -- empty list =#,
        [] #=run -- empty list =#,
        [] #=isrunning=#)
end

duffing

Test to see if it generates a window with a single “ball” shown.

d1=duffing()

UndefVarError: `ax` not defined

Stacktrace:
 [1] duffing(; xv::Vector{Float64}, p::dp, t::Float64, dt::Float64, delay::Float64, step::Function, dxdt::Function)
   @ Main ./In[4]:14
 [2] duffing()
   @ Main ./In[4]:1
 [3] top-level scope
   @ In[5]:1

Next, see if you can move the ball around with

d1.ball[]=Point2f(0.2, 0.3)

## Activity 2

We now need a method which updates xv and the ball object. This is analogous to the setangle! method from the Animating a Pendulum lecture. Test it.

function setxv!(d::duffing,xv)
    # fill in details
    #
end

## Activity 3

Now we want to modify our object so that we can click on the image to update xv. Copy you constructor from above, and add the necessary code (using the select_point function) so that it listens for a mouse click, and does the appropriate updating. Look in the Animating the Pendulum notes.

Don’t forget to disable the rectanglezoom function. Now might be a good time to add code that sets the x and y axis limits.

Note: you want to use the setxv! command that you just wrote. Thus you need to set up the interactive elements after you instantiate the object

## Activity 4

Now lets try to implement the timestepping. Make a method step! which will use the integrator to take one timestep. It should call setxv! – it should also update the time. This is important, because we have a time dependent force.

function step!(d::duffing)
    # fill in code
    # Careful -- look at the arguments of d.step
    # -- they are different from those in the pendulum case
end

Debug, and see if it works for single steps.

step!(d1)  # or whatever else you called your instance

Then try a whole sequence of steps. Click on a few points, and then run this loop a few times. The trajectories should make sense.

or i in 1:1000
    step!(d1)
    sleep(d1.delay)
end

## Activity 5

Now we want to add the run button. Copy your constructor code to below here. Add the code necessary to add a run button to the duffing object. The key lines are:

isrunning = Observable(false)
label = map(cond -> cond ? "Stop" : "Run", isrunning)
run = Button(fig[2,1]; label = label, tellwidth = false)
on(run.clicks) do clicks; isrunning[] = !isrunning[]; end
startsim(d)

You will also need the startsim method

function startsim(p::duffing)
    isrunning_notifier = Condition()
    on(cond -> cond && notify(isrunning_notifier), p.isrunning)
        errormonitor(@async while true
        if p.isrunning[]
            isopen(p.fig.scene) || break # stops if window is closed
            step!(p)
            sleep(p.delay) 
        else
            wait(isrunning_notifier)
        end
    end)
end

Note: debugging this can be hard, as the code inside the @async block does not throw very good error messages. (Or any at all. I think wrapping it in errormonitor, as I did, should turn the error messages back on, but I haven’t played with it enough to be sure.) If it doesn’t work, your best bet is to ask for some help with debugging, as it is hard to give generic advice here. Make double-sure that step! works properly – as that is the most likely place where there is a bug.

Add this code to your constructor, and store run and isrunning in your object.

Test it out. Click on different points, and watch the trajectory evolve. Do the results make sense? (Keep F=0 so that it is easy to interpret)

## Activity 6

We now want to add sliders to be able to adjust F and ω. Here is an example of how to add sliders to an existing duffing object:

d3=duffing()
sg=SliderGrid(d3.fig[1,2],
    (label="F",range=0.:0.01:1.,format = "{:.2f}",startvalue=0.),
    (label="ω",range=0.01:0.01:2,format = "{:.2f}",startvalue=1.),
    width=200,tellheight=false)

Those sliders don’t do anything yet. Feel free to drag them around, but it will not make any change to the simulation.

We then need to make it so that when the sliders are adjusted, the attributes of the oscillator are updated

Fobserver=sg.sliders[1].value;
ωobserver=sg.sliders[2].value;

on(Fobserver) do Fobserver
    d3.p.F=Fobserver[]
end

on(ωobserver) do ωobserver
    d3.p.ω=ωobserver[]
end

Play with the F slider. If it runs too slow or too fast, you can adjust the speed by changing the timestep “on the fly” with something like the following line

d3.dt=0.04

Again, copy your constructor, and edit it so that it automatically adds the sliders – and stores them in the object. Test. (Note: You don’t need the d3.p or d3.fig. Inside the constructor you have direct access to p or fig. You do not want d3 to appear anywhere – as that is a particular instance of the object, not the one we are creating.)

Make sure to test your code.

## Activity 7

Lets now add the trails. For this you will need to edit both the constructor, and the step! method. Look up how we did this in Animating a Pendulum, and implement it.

## Activity 8

Copy your code here, so that in one cell you have the whole stand-alone program. Restart the kernel and run the cell to make sure it work. That way if you want to show anyone your simulation, you can just go to this one cell, and not have to run a bunch of different things.

Play with your simulator.

Find parameters where you have regular (non-chaotic) trajectories. Watch the behavior.

Find chaotic parameters. Watch the behavior.

We will go over this a bit in class, and try to understand what is happening.

## Activity 9: Bonus

This is a bonus activity, that you don’t need to do, but is fun.

Modify your code so that it tracks two oscillators. Make it so that when you click on the graph, one of them gets set to that point, and the other gets set to a slightly offset point. That way you can observe how close the trajectories stay together. Make the points and trails different colors.

An even better way to get insite is to not just 2, but \(n\) oscillators, where \(n\) can be set at instantiation time. When you click on a point, these get initiated to a ring of points around where you click. Take \(n~1000\) to get an understanding of how chaos works in this system