The Duffing Oscillator is a mass$-$spring system with non$-$linear restoring force

$F_{\text{r}\text{e}\text{s}\text{t}\text{o}\text{r}\text{i}\text{n}\text{g}}=-{}_0^{2}y-y^3.$ Assuming a mass of $1$ kilogram, its motion can be described by

a second order ODE derived from Newton's second law of motion. Taking into account

damping, and an oscillatory driving term, the equation of motion is

$y^{}+y^{}+{}_0^{2}y+y^3=f_{0}sin(wt),$

where $y$ is the distance from equilibrium, and $,{}_0,,f_0$ and $w$ are treated as parameters

of the model.

The ODE is nonlinear and thus difficult to solve analytically. Numerical techniques

however, are well suited to analysing the system.

a$)$ Setting $v_1=y,v_2=y^{},$ rewrite the ODE $(2)$ as a system of first order ODEs of the

form:

$v_1^{}=f(t,v_1,v_2),$

$v_2^{}=g(t,v_1,v_2).$

Copy the following code into MATLAB, and name it osc.m$.$

Modify lines $26$ and $27$ to reflect the formulas you've derived in equations $(3)$ and

$(4)$ above. Fix the parameter values ${}_0=1,=0\mathrm{.}2,=0\mathrm{.}03$ and $w=1\mathrm{.}2$ in lines

$13$ through $16.$

Submit the modified version of osc.m as part of your assignment.

function $f=$osc$(t,y,F0)$

$\%$ Return a column vector

$f=[0$;$0]$;

$\%$ Oscillator Parameters

$\%$ omega$0-$ resonant frequency

$\%$ gamma $-$ damping

$\%$ eta $-$ nonlinearity, set to $0$ for harmonic

$\%$ w $-$ drive frequency

$\%$ FO $-$ drive force

omega $0=8$ Paste value here;

gamma $=$ opaste value here;

eta $=\%$ Paste value here;

$w=$ oraste value here;

$\%v1=y,v2=$ydot

$v1=y(1)$;

$v2=y(2)$;

$\%$ Force Calculation

$F=F0**sin(w****t)$;

$\%$system of First Order ODEs

$f(1)=$ opaste formula for $v_{-1}$ dot in here

$f(2)=$ opaste formula for $v_{-}2$ dot in here.

end $($b$)$ Let$$s consider the following two driving amplitudes:

$($i$)$ f$0=0\mathrm{.}4$

$($ii$)$ f$0=0\mathrm{.}7.$

Create an m$-$file to produce plots of solution curves to the two driving amplitudes.

$$ Your code should call osc.m $($found in part a$)),$ and ode$45.$

$$ Set initial conditions to be y$(0)=$y$(0)=0.$

$$ The domain of t should be $[0,1000].$

$$ Name the m$-$file you used plotsoln.m$,$ and submit it as part of your assignment.

In both cases, your solution should be oscillatory and should $$settle down into its

final rhythm$$ after about t $600$ or so$.$

What is the approximate amplitude of the solution$$s final waveform for the two

driving amplitudes $($we$$ll call this the $$final amplitude$)?$

For each value of f$0,$ submit plots of the solution curves.

$($c$)$ We now investigate the relationship between f$0$ and the waveform$$s final amplitude.

Write another program called drive.m$.$ This code should find the waveform$$s final

amplitude for $20$ different values of f$0$ between $0\mathrm{.}4$ and $0\mathrm{.}7.$

Submit a plot which summarises this data. How does the final amplitude respond

to changes in f$0?$

Submit the plot, and drive.m as part of your assignment.