Parametric resonance

(Go down and play with controls; Requires Mathematica 6 )

We solve the equation of motion

x''[t]+2λ x'[t]+((ω0)^2)(1+h Cos[ω t])x[t]==0

with ω0=1 in the time interval 0 ≤ t ≤ tMax. The instability increment is given by

s=(1/2)((Re[(((((h ω0)/2)))^2)-(((ω0-2ω))^2)])^(1/2))-λ

and it reaches a maximum at ω = ω0/2. For the instability to occur, pumping h should prevail over damping, s>0.

Plotting phase trajectories

PlotPhaseTrajectory[ω0_,ω_,h_,λ_,x0v0_,tMax_]:= (* This subroutine solves the equation of motion and plots trajectories, depending on input parameters *)
Quiet[Module[{nds,x},

Eq=x''[t]+2λ x'[t]+((ω0)^2)(1+h Cos[ω t])x[t]==0;
ω0=1;
InstabilityIncrement=(1/2)((Re[(((((h ω0)/2)))^2)-(((ω0-2ω))^2)])^(1/2))-λ;
(* Print["InstabilityIncrement = ",InstabilityIncrement]; *)

nds =NDSolve[Join[{Eq},Thread[{x[0],x'[0]}==x0v0]],{x},{t,0,tMax},PrecisionGoal->6];

(* Now plotting the phase trajectory *)
ParametricPlot[{x[t],x'[t]}/.nds[[1]],{t,0,tMax}, 
Frame->True,Axes -> False, Prolog -> {GrayLevel[.5], PointSize[0.02], Point[{0,0}]},ImageSize->{600,600},PlotRange->{{-3,3},{-3,3}},
FrameTicks -> None,
PlotPoints -> Ceiling[8 tMax+4],MaxRecursion->ControlActive[2, Automatic],
PlotStyle->{Thickness[.005]}] 
]]

Manipulate[
PlotPhaseTrajectory[ω0,ω,h,λ,x0v0,tMax],
{{tMax,50,"maximal time"}, 0, 200},
{{ω,2ω0,"frequency of perturbation ω"},0,4, Appearance -> "Labeled"},
{{h,0.1,"amplitude of perturbation h"},0,0.5, Appearance -> "Labeled"},
{{λ, 0,"damping λ"},0,0.3, Appearance -> "Labeled"},
{{x0v0,{1, 1}, "initial\nposition\n{x0,v0}"},{-2,-2},{2,2},ImageSize->Medium},
ControlPlacement->{Top,Top,Top,Top,Left},SaveDefinitions->True,AutorunSequencing->{4,5},TrackedSymbols:>{tMax,ω,h,λ,x0v0}]

Plotting the energy of the oscillator

Energy[ω0_,ω_,h_,λ_,x0v0_,tMax_]:= (* This subroutine solves the equation of motion and plots trajectories, depending on input parameters *)
Quiet[Module[{nds,x},

Eq=x''[t]+2λ x'[t]+((ω0)^2)(1+h Cos[ω t])x[t]==0;
ω0=1;

nds =NDSolve[Join[{Eq},Thread[{x[0],x'[0]}==x0v0]],{x},{t,0,tMax},PrecisionGoal->6];

e[t_]:=(((x'[t])^2)+((ω0)^2)((x[t])^2))/.nds[[1]];   (* Proportional to the energy of the oscillator *)
Plot[e[t],{t,0,tMax}]

]]

Manipulate[
Energy[ω0,ω,h,λ,x0v0,tMax],
{{tMax,50,"maximal time"}, 0, 200},
{{ω,2ω0,"frequency of perturbation ω"},0,4, Appearance -> "Labeled"},
{{h,0.1,"amplitude of perturbation h"},0,0.5, Appearance -> "Labeled"},
{{λ, 0,"damping λ"},0,0.3, Appearance -> "Labeled"},
{{x0v0,{1, 1}, "initial\nposition\n{x0,v0}"},{-2,-2},{2,2},ImageSize->Medium},
ControlPlacement->{Top,Top,Top,Top,Left},SaveDefinitions->True,AutorunSequencing->{4,5},TrackedSymbols:>{tMax,ω,h,λ,x0v0}]