Math 3280 - Hampton - Lab 6- Due 7/17

Resonance

In this part of the lab we will consider problem number 27 from Chapter 5.6 in the text.  For convenience, we will use m= k = 1.

UnitMassUnitSpringODESolution[c_,w_]:=DSolve[{y''[t]+c y'[t]+y[t]==Cos[w t]},y[t],t][[1]]

RandomParameters={C[1]->Random[],C[2]->Random[],c->1+Random[]*3,w->1+Random[]*9}

ExampleSolution[t_]:=Simplify[y[t]/.UnitMassUnitSpringODESolution[c,w]/.RandomParameters]

We are only interested in the steady state behavior of the forced oscillator.  Here is an example of the complete solution decaying to the steady-state.

ExampleSolution[t]

The DSolve command sometimes generates complex answers with a very small imaginary part; these are errors.  In the plot command below we extract the real part of the solution using the Re[] command.

Plot[Evaluate[Re[ExampleSolution[t]]],{t,0,30}]

Lets extract the steady-state solution.

Simplify[y[t]/.UnitMassUnitSpringODESolution[c,w]/.{C[1]->0,C[2]->0}]

The total amplitude is the square root of the sum of the squares of each amplitude:

Simplify[((w^2-1)^2+(c w)^2)^(1/2)/(1+(-2+(c^2)) (w^2)+(w^4))]

So finally we define the SteadyStateAmplitude function. This corresponds to C(w) in the text.

SteadyStateAmplitude[c_,w_]:=(1/(Sqrt[1+(-2+(c^2)) (w^2)+(w^4)]))

Exercises

1. Plot the SteadyStateAmplitude function in the c-w plane using the functions Plot3D and ContourPlot, on the domain c in [0,4] and w in [0,8].  To get an accurate picture, use the ViewPoint, PlotPoints, and PlotRange options.

2. Make a plot of the SteadyStateAmplitude versus c along the curve of maximum amplitude (i.e., for each value of c, choose the value of w which maximizes the SteadyStateAmplitude function) for c in [0,4].  Be careful - the curve of maximum amplitude is continuous, but not always differentiable in this range of c values.

Lunar mission example - based on example from pages 468 - 469 in the text.

First read the example on page 468-469.  Then do the exercises below using the function DrawOrbit below.  The arguments of DrawOrbit are the initial position and velocity in the x1 and x2 directions, and the ending time of the numerical approximation.  You may want to alter some of the settings in the function or define your own variant of it.

Mu=12277471/1000000000;  (* Mu is entered as an exact fraction to encourage mathematica to use more precision *)
DrawOrbit[x10_,x1v0_,x20_,x2v0_,time_]:=Module[
{sol,ptfun,pttab (* temporary variables for this module *)},sol={x1,x2}/.NDSolve[{x1''[t]==x1[t]+2x2'[t]- (1-Mu)(x1[t]+Mu)/((x1[t]+Mu)^2+(x2[t])^2)^(3/2)-(Mu)(x1[t]-1+Mu)/((x1[t]-1+Mu)^2+(x2[t])^2)^(3/2),x2''[t]==x2[t]-2x1'[t]-(1-Mu)(x2[t])/((x1[t]+Mu)^2+(x2[t])^2)^(3/2)-(Mu)(x2[t])/((x1[t]-1+Mu)^2+(x2[t])^2)^(3/2),x1[0]==x10,x1'[0]==x1v0,x2[0]==x20,x2'[0]==x2v0},{x1,x2},{t,0,time},AccuracyGoal->18,MaxStepSize->.005,MaxSteps->Infinity][[1]];
ptfun[ti_]:={Evaluate[sol[[1]][ti]],Evaluate[sol[[2]][ti]]};
pttab=Table[ptfun[t],{t,0,time,Min[time/200,.5]}];
Show[Graphics[{Point[{-Mu,0}],Point[{1-Mu,0}],Line[pttab]}],PlotRange->All]
]

DrawOrbit[.8,-.01,0,.01,20]

Exercises

3. Re-create the example pictures on page 469 using the DrawOrbit function or a variant of it.

4.  Find an interesting orbit that stays near the point x1 = 1/2 - Mu, x2 = (Sqrt[3])/2 for a while and then comes close to both the Earth and the Moon.