Pendulum

How to integrate with a computer

AP E3200x

Introduction

The equation of motion for a simple pendulum is a good way to illustrate how a computer can solve a differential equation.

The position of a pendulum is given by the angle θ. Then,
	(∂_t) θ = ω(t)
	(∂_t) ω = - (g/L) Sin[θ]
when the angle is small |θ| << 1, then Sin[θ] ~ θ, and
	(∂_t) (∂_t) θ = - (g/L) θ
θ(t) = (Subscript[θ, 0])Sin[(Subscript[ω, 0])t + (Subscript[ϕ, 0])] where ((Subscript[ω, 0])^2) = g/L, and (Subscript[θ, 0]) and (Subscript[ϕ, 0]) are initial conditions.

Integration with a Computer

The numerical solution for the pendulum's motion requires the integration of two coupled first order ordinary differential equations (ODEs). This is the best problem to begin the exploration of computation physics. 

We will first explore Euler's method: integration of the equations of motion by approximating the first-derivative with a forward difference. We'll see that Euler's method fails and becomes unstable. 

Next, we'll explore the Euler-Cromer method which significantly improves the stability of the integration. It is similar to the Leap-Frog method which is based on central-difference approximation to the first-derivative.

Finally, we'll look into the numerical tools built-in to Mathematica.

dt = 0.01; MachineNumberQ[dt]

Euler's Method

Using two forward-difference approximations…

stepE[{θ_, ω_}] := N[{θ + ω dt, ω - Sin[θ] dt}]

stepE[{θ0, ω0}]

stepE[{N[75 Degree],0.0}]

(steps = NestList[stepE, {N[75 Degree],0.0}, 5000];) // Timing

ListPlot[steps]

ListPlot[First[Transpose[steps]],PlotLabel->"θ vs. time step"]

If dt gets too large, then the entire method fails.

dt = 0.1;steps = NestList[stepE, {N[75 Degree],0.0}, 5000];

ListPlot[First[Transpose[steps]],PlotLabel->"θ vs. time step"]

Central Difference (Leap-Frog Method)

The leap-frog method approximates the differential equations with central-difference approximations.
	(∂_t) θ ~ ((θ[t+δt/2]  -  θ[t-δt/2])/(δt)) = ω(t) + O(((δt)^2))
	(∂_t) ω ~ ((ω[t+δt/2]  -  ω[t-δt/2])/(δt))= -  Sin[θ(t)]+ O(((δt)^2))
Looking carefully at the difference equations, evaluation of θ steps "leap" over ω, and evaluations of ω steps "leap" over θ. This presents two problems: first, the method is not self-starting, and, second, the values of θ and ω are not known at the same time step. This method has the advantage of being second-order with reduced truncation error.
In the following, we'll use an Euler "half-step" to start the integration by computing θ(δt/2). Then, we'll "leap" with ω followed by θ using "full steps".

stepLF[{θ_, ω_}] := Module[{ωnext},
ωnext = N[ω - Sin[θ] dt];
{N[θ + dt ωnext], ωnext}]

dt = 0.01;
w0 = 0.0;
θ0 = N[75 Degree];
θ1 = N[θ0 + w0 (dt/2)]; (* Euler Step to start *)(steps = NestList[stepLF,{θ1, w0}, 5000];) // Timing

ListPlot[First[Transpose[steps]],PlotLabel->"θ vs. time step"]

Euler-Cromer Method

This method is practically the same as the leap-frog method. The only difference is to drop the pretense of leaping and shift the two variables, {θ, ω}, so that they "align" (instead of "stagger") every time step. In Euler-Cromer, we actually make "two steps": first a step in the frequency, ω, and then in the angle, θ; however, we do not need to "start" the integration with a "half-step".

stepEC[{θ_, ω_}] := Module[{ωnext},
ωnext = N[ω - Sin[θ] dt];
{N[θ + dt ωnext], ωnext}]

dt = 0.01;(steps = NestList[stepEC, {N[75 Degree],0.0}, 5000];) // Timing

(Mathematica requires more time when using the Euler-Cromer method, but this can be fixed by pre-compiling the module. We'll discuss this later in the course.)

ListPlot[steps]

ListPlot[First[Transpose[steps]],PlotLabel->"θ vs. time step"]

The Euler-Cramer method is also more stable for large time steps.

dt = 0.1;steps = NestList[stepEC, {N[75 Degree],0.0}, 5000]; // Timing

ListPlot[First[Transpose[steps]],PlotLabel->"θ vs. time step"]

ListPlot[First[Transpose[steps[[Range[1,500]]]]],PlotLabel->"θ vs. time step"]

Using NDSolve

Mathematica can numerically integrate differential equations using NDSolve[...]. NDSolve is probably the most sophisticated and comprehensive command within Mathematica. Since there are many methods to finding numerical solutions to differential equations, NDSolve has many options.

In order to find the solution to the equations of motion for the pendulum, we need to define:
	1. The unknown functions of time, θ[t] and ω[t].
	2. The equations of motions, and
	3. The equations for the initial conditions.

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

initial = {θ[0] ==N[75 Degree], w[0] == 0}

?NDSolve

Options[NDSolve]

(sol = NDSolve[ eqs ~ Join ~ initial, {θ[t], w[t]}, {t, 0, 60}] )// Timing

Plot[θ[t]/.sol,{t,0,60}]

To appreciate the numerical solution, examine

θ[t] /. sol// FullForm // Shallow

The interpolating function contains all of the numerical solution. For example, the time points are found by

Part[ θ[t] /. sol , 1, 0, 3, 1] // Shallow

stepsND = (Part[ θ[t] /. sol , 1, 0, 3, 1] // Length)

ListPlot[Delete[RotateLeft[(θ[t]/.sol)[[1,0,3,1]]]-(θ[t]/.sol)[[1,0,3,1]],-1],PlotRange->All,PlotLabel->"Time Steps for NDSolve"]

Each solution can be viewed continuously by plotting the smooth curve between the points.

Other Integration Methods in NDSolve

(solRK = NDSolve[ eqs ~ Join ~ initial, {θ[t], w[t]}, {t, 0, 60}, Method->ExplicitRungeKutta] )// Timing

Plot[Evaluate[θ[t]/.sol]-Evaluate[θ[t]/.solRK],{t,0,60}]

stepsND = (Part[ θ[t] /. solRK , 1, 0, 3, 1] // Length)

ListPlot[Delete[RotateLeft[(θ[t]/.solRK)[[1,0,3,1]]]-(θ[t]/.solRK)[[1,0,3,1]],-1],PlotRange->All,PlotLabel->"Time Steps for Runge-Kutta"]

Understanding Interpolating Function

As an example of using interpolating functions, consider the following.

sineList = Table[{t, Sin[t]}, {t, 0, 10, 0.4}];

Length[sineList]

lp=ListPlot[sineList,PlotStyle->{PointSize[0.02`]}]

sineSmooth = Interpolation[sineList]

sp=Plot[sineSmooth[t],{t,0,10}]

Show[lp,sp]

Show[lp,sp,PlotRange->{{0,2},All}]

Interpolating functions are "real" functions that can be differentiated and integrated.

In NDSolve[...], the number of points determines the how long the numerical solution can be solved with a given accuracy. The default value is 1000 points, and for some complex problems we will need to increase this number by at least a factor of 10!

Comparing the NDSolve with Euler-Cromer's Method

dt = 0.1;steps = NestList[stepEC, {N[75 Degree],0.0}, 500];

solEC = Table[{(i-1) dt, First[steps[[i]]]}, {i, 1, 500}];

pEC=ListPlot[solEC,PlotStyle->{RGBColor[0,0,1]}]

pNDSolve=Plot[θ[t]/.sol,{t,0,500 dt}]

Show[pEC,pNDSolve]

Looking at the difference...

solND[tt_] = θ[t] /. First[sol]/. t-> tt;err = Table[Flatten[{(i-1) dt, First[steps[[i]]] - solND[(i-1) dt]}], {i, 1, 500}];

ListPlot[err]

Nonlinear Pendulum

Finally, we can calculate how the pendulum’s period changes with amplitude.

initial = {θ[0] ==N[42 Degree], w[0] == 0}
sol = NDSolve[ eqs ~ Join ~ initial, {θ[t], w[t]}, {t, 0, 60}];

Plot[θ[t]/.sol,{t,0,10}]

FindRoot[θ[t]/.sol,{t,2}]

tRoot = t /. %

findFrequency[initialθ_] := Module[{sol, initial,tNew},
initial = {θ[0] ==N[initialθ Degree], w[0] == 0};
sol = NDSolve[ eqs ~ Join ~ initial, {θ[t], w[t]}, {t, 0, 100}];
tNew = t /.  FindRoot[θ[t]/.sol,{t,tRoot}];
tRoot = tNew;
π/(4tNew) ]

findFrequency[1]
findFrequency[75]
findFrequency[175]

ListPlot[Table[{angle, findFrequency[angle]}, {angle, 1, 179, 0.5}], Frame -> True, Joined -> True,FrameLabel -> {"Initial θ[0] (Degree)","Pendulum Frequency"}, PlotLabel -> "Nonlinear Pendulum"]

initial = {θ[0] ==N[179 Degree], w[0] == 0};
sol179= NDSolve[ eqs ~ Join ~ initial, {θ[t], w[t]}, {t, 0, 80}];initial = {θ[0] ==N[155 Degree], w[0] == 0};
sol155 = NDSolve[ eqs ~ Join ~ initial, {θ[t], w[t]}, {t, 0, 80}];
initial = {θ[0] ==N[105 Degree], w[0] == 0};
sol105 = NDSolve[ eqs ~ Join ~ initial, {θ[t], w[t]}, {t, 0, 80}];
initial = {θ[0] ==N[75 Degree], w[0] == 0};
sol75 = NDSolve[ eqs ~ Join ~ initial, {θ[t], w[t]}, {t, 0, 80}];
ParametricPlot[{{θ[t],w[t]}/. sol75,{θ[t],w[t]}/. sol105,{θ[t],w[t]}/. sol155,{θ[t],w[t]}/. sol179},{t,0,80}, AxesLabel -> {"θ[t]","w[t]"}, PlotLabel -> "Nonlinear Pendulum"]

The nonlinear pendulum can not be described by Sine and Cosine functions. (But, they can be expressed in terms of Elliptic functions.)

Summary

Computing the solution to ordinary differential equations is the most fundamental activity of computational physics. This is because Newton's Laws of motion are expressed as differential equations. 

A computer can only approximate the actual solution. A computer finds a list of coordinates that define a trajectory. Typically, we start with an initial position and velocity, and compute the next position and velocity using a finite-difference approximation to the time derivatives. This is not a trivial task. 

Euler's method fails for reasonably long integration times. 

Euler-Cromer's method is both more accurate and more stable. 

NDSolve is built into to Mathematica. NDSolve uses one of several advanced schemes with "adaptive" step sizes to achieve a specific accuracy. NDSolve is also relative fast, and, for most of our physics experiments, NDSolve is the method of choice for both accuracy and speed.