Assignment 5

Name:
Student Number:

Q1

Remove["Global`*"]

Numerical solution of the Duffing oscillator.  See onilne notes 7.5 and 7.9

The Duffing oscillator is described by the equation

(((d^2)x)/(dt^2))+γ(dx/dt)+(Subscript[ω, 0]^2)x+((βx)^3)=F		(Duffing oscillator)

where F is oscillatory forcing term.  We will use F(t) = F0 sin(ωf t).
Solve numerically for x(t) with the following parameters

xi=1.0; (* initial displacement *)
 ui=0; (* initial velocity *)

γ=0.01; (* damping frequency*)
ω0=1.0; (* natural frequency of linear, undamped version of oscillator *)
β=1; (* strength of the nonlinear term *)

ωf=ω0/4;  (* driving frequency *)
F0=5.0; (* strenght of driving term *)

T=200;  (* Simulation period *)
Δt=0.1; (* Time step *)

using Euler, 2nd order RK and 4th order RK.  For each method, plot the resulting interpolated solution x(t), as well as u (velocity) as a function of x parametrically.

a) Euler - Euler will fail with Δt=0.1.  Try again with smaller Δt until you get something stable

b) 2nd order RK

c) 4th order RK

Q2

Solve the same Duffing oscillator probem as in Q1, but implement the following predictor-corrector algorithm

ytmp = (Subscript[y, n]) + h f((Subscript[y, n]), (Subscript[t, n]))
(Subscript[y, n+1]) = (Subscript[y, n]) + (1/2)h ( f((Subscript[y, n]), (Subscript[t, n])) + f(ytmp, (Subscript[t, n+1])) )

which can also be written

s1=f((Subscript[y, n]), (Subscript[t, n]))
ytmp = (Subscript[y, n]) + h s1
s2=f(ytmp, (Subscript[t, n+1]))
(Subscript[y, n+1]) = (Subscript[y, n]) + (1/2)h ( s1 + s2 )

Δt = 0.1 should get you a very similar result to 2nd order RK.  Decrease Δt until you get a similar result to RK4.

Q3

The famous Lorenz equations represent a simplified model of convection in the atmosphere.  They are three coupled, non-linear first order ODE’s for x(t), y(t) and z(t):

(dx/dt)=σ ( y-x)
(dy/dt)=x (ρ -z)-y
(dz/dt)=x y - β z

where σ, β and ρ.  Use RK4 to solve these equations with

σ=10;
β=8/3;
ρ=28;

xi=1;
yi=0;
zi=0;

T=50;
Δt=0.01;

You will need to adapt the algorithm to solve three simulataneous ODE’s.

Plot your interpolated solutions parametrically using ParametricPlot3D, for which you will need to use the option PlotPoints->1000 (or similar) to get a nice smooth plot. Combine your plot of the trajectory with ListPointPlot3D and Manipulate to get an animation of the (x, y, z) point moving in time along the trajectory.  In using Manipulate, make sure your increment in time is small enough, e.g., by using {t, 0, T, Δt/2} in specifying the time values over which to animate.