Weekly Exercises (2)

Assigned: January 29, 2013
Due: before class on February 7, 2013

Verlet's Method (in-class)

Verlet's method uses a central-difference for the second-derivative of θ:
	(∂_t) (∂_t) θ = - (g/L) Sin[θ] , or approximately as:
	((θ[t-δt] - 2θ[t] + θ[t+δt])/(((δt))^2))~ -Sin[θ[t]] + O(((δt)^2))
This method computes only the position and not the velocity. Consequently, it is suitable for problems when the force depends only on the position (and not the velocity). Notice that the trunction error scales like ((δt)^4). Verlet's method is relatively accurate and easy-to-program. Like the leap-frog method, this method is not self-starting. Use one step of Euler to start your integration, and then compute an approximation to the angular motion of the pendulum, θ(t), for several periods. You may start with any initial condition. You may also normalize length and time so that (g/L) = 1.0.

A Solution (Method 1)

stepV[θ_]:= Module[{tlast},
tlast = thetaLast;
thetaLast = θ;
-dt^2 Sin[θ] + 2 θ - tlast]

dt = 0.01;
w0 = 0.0;
θ0 = N[75 Degree];
thetaLast = θ0 - w0 dt;
steps = NestList[stepV, θ0, 5000];

steps // Shallow

ListPlot[steps]

A Solution (Method 2)

stepV[{θNow_,θThen_}]:= {-dt^2 Sin[θNow] + 2 θNow - θThen, θNow}

dt = 0.01;
w0 = 0.0;
θ0 = N[75 Degree];
thetaLast = θ0 - w0 dt;
steps = NestList[stepV, {θ0 ,thetaLast}, 5000];

steps // Shallow

ListPlot[First[Transpose[steps]]]

Period of a Nonlinear Pendulum (in-class)

A linear pendulum has a period that is independent  of the amplitude of oscillation. This is why pendulums can be used for clocks. A real pendulum is nonlinear, and the period changes as the amplitude of its swinging increases.

Use NDSolve[…] and compute the period of our pendulum as a function of angle.
You should begin with the solution discussed in class:

eqs = {θ'[t] == w[t], w'[t] == - Sin[θ[t]]}

initial={θ[0]==N[75 °],w[0]==0};
sol=NDSolve[Join[eqs,initial],{θ[t],w[t]},{t,0,20}];
Plot[Evaluate[θ[t]/.sol],{t,0,20}]

Next, you will need to use FindRoot[…], to estimate the time when the pendulum crosses the horizontal. This will be 1/4 of the period. FindRoot[...] returns a "Replacement Rule" (not a number). We must use "/." (called "slash-dot", and meaning "ReplaceAll") to apply the rule.

?FindRoot

FindRoot[Evaluate[θ[t] /. sol], {t, 2.0}]

(θ[t] /. sol)/. % (* should be near zero *)

Note: FindRoot[...] needs a "guess" where to look for a zero-crossing. The period of the pendulum is four times the time to the first "zero crossing".

4 t /. FindRoot[Evaluate[θ[t] /. sol], {t, 2.0}]

Combine these statements into a single Module[…] and define a function of the initial angle that returns the period.

A Solution

period[θ0_] := Module[{initial,sol},
initial = {θ[0] ==θ0, w[0] == 0};
sol = NDSolve[ eqs ~ Join ~ initial, {θ[t], w[t]}, {t, 0, 20}] ;
4 t /. FindRoot[Evaluate[θ[t] /. sol], {t,2.0}]]

pList = Table[{t0/Degree, period[t0]/(2 π)}, {t0, 1 Degree, 170 Degree, 5 Degree}];

ListPlot[pList]

A "Constant" Restoring Force

Compute the motion of an object confined by a "constant" restoring force. The force is directed towards the origin, x = 0. When x < 0, then the force is (Subscript[F, 0]); and when x > 0, then the force is – (Subscript[F, 0]). (For simplicity, let (Subscript[F, 0]) = 1.0.)
Use both the Euler-Cromer method and the built-in method in Mathematica, NDSolve[…]. Start with the object at rest (v = 0.0) but with a non-zero initial position, x(0) ≠ 0. Note how the period of oscillation depends upon the amplitude. Finally, use NDSolve[…] and the technique used in the previous exercise to graph the period of oscillation as a function of the initial position. (Did you notice the geometric shape of the trajectory, x(t)? The period should scale as (Sqrt[x(0)]).)

Solution

First, the Euler-Cromer method:

stepEC[{x_, v_}] := Module[{vnext},
vnext = If[x < 0.0,v +  1.0 dt, v -1.0 dt];
{N[x + dt vnext], vnext}]

dt=0.01;
steps=NestList[stepEC,{1.,0.},5000];
ListPlot[First[Transpose[steps]],PlotLabel->"x vs. time step"]

Then, we use NDSolve[…]:

eqs = {x'[t] == v[t], v'[t] == If[x[t] < 0, 1.0, -1.0]}

initial={x[0]==1.,v[0]==0.};
sol=NDSolve[Join[eqs,initial],{x[t],v[t]},{t,0,30}];
Plot[Evaluate[x[t]/.sol],{t,0,30}, AxesLabel -> {"t","x[t]"}]

initial={x[0]==5.,v[0]==0.};
sol=NDSolve[Join[eqs,initial],{x[t],v[t]},{t,0,30}];
Plot[Evaluate[x[t]/.sol],{t,0,30}, AxesLabel -> {"t","x[t]"}]

The period gets longer as the amplitude gets larger. We also notice that the trajectories are parabolas!

Finally, we compute the period as a function of initial position:

period[x0_] := Module[{initial,sol},
initial = {x[0] == x0, v[0] == 0.0};
sol = NDSolve[ eqs ~ Join ~ initial, {x[t], v[t]}, {t, 0, 40}] ;
4*(t /. FindRoot[Evaluate[x[t] /. sol], {t, 2.0}])]

pList = Table[{x0, period[x0]},{x0, 1.0,5.0, 0.1}];

p1=ListPlot[pList,PlotLabel->"Period vs. Initial Position",PlotStyle->PointSize[0.02], AxesLabel -> {"(Subscript[x, 0])","T"}]

And we show that the period scales like (Sqrt[x(0)]),

Fit[pList, Sqrt[x], x]

p2=Plot[5.65685 (Sqrt[x]),{x,1.,5.},PlotStyle->{RGBColor[0,0,1]}];
Show[p1,p2]

Do you know the exact relationship between oscillation period and (Sqrt[(Subscript[x, 0])])?