Spring 1998

Maple Lab, Week 26

The Pendulum Equation and Maple’s dsolve Command

Background: The second order differential equation

(1)  \frac{d^{2}θ}{dt^{2}} + c\frac{dθ}{dt} + \frac{g}{l} sin(θ) = 0

is a mathematical model for the motion of a pendulum. The pendulum consists of a rigid, weightless rod of length l, pinned at one end, and with a mass m at the other end. The rod pivots about the pinned end, O, and is assumed to move in a single plane. The angular position of the rod at time t, measured counterclockwise from vertically downward, is θ = θ(t). The term c\frac{dθ}{dt} accounts for the drag on the mass due to air resistance. You will be asked to derive the equation (1) in the exercises. If we simplify (1) by assuming \frac{g}{l} = 1 and c = 0, we get

(2)  \frac{d^{2}θ}{dt^{2}} + sin(θ) = 0.

Maple has a command, dsolve, for solving differential equations. The following example illustrates the use of dsolve to solve the simple initial value problem y′ +3y = 0, y(0) = 2. Additional examples may be found in the online help.

> de1 := diff(y(t),t) + 3*y(t) = 0;
> init1 := y(0) = 2;
> sol1 := dsolve({de1,init1},y(t));

Note that the solution is an equation. If we would like to plot the solution, for example, we could use the plot command together with the rhs command (for right hand side):

> ysol1 := rhs(sol1);
> plot(ysol1, t=0..1);

Sometimes dsolve can’t solve a differential equation. For example, the equation (2) is nonlinear and happens to be too difficult for Maple, as the following commands illustrate:

> de2 := diff(th(t),t,t) + sin(th(t)) = 0;
> init2 := th(0) = 0, D(th)(0) = 2;
> sol2 := dsolve({de2,init2},th(t));

There are at least two things which can be done to remedy the situation: simplify the differential equation, or seek an approximate (i.e., numerical) solution. As an example of the first approach, we can approximate sin(θ) by the linear approximation L(θ) = θ near θ = 0. This gives the equation θ′′ + θ = 0, which is easy for Maple to solve (and easy to solve by hand, as well):

> de3 := diff(th(t),t,t) + th(t) = 0;
> sol3 := dsolve({de3,init2},th(t));

To obtain a numerical solution, include type = numeric, or simply numeric, in the calling sequence:

> sol2 := dsolve({de2,init2},th(t),numeric);

The solution is a procedure which can be evaluated at any value of the independent variable:

> sol2(0); sol2(1);

One way to plot the solution is to create a list of ordered pairs and use the plot command:

> plt := NULL:
> for i from 0 by .1 to 10 do
> plt := plt, [rhs(sol2(i)[1]), rhs(sol2(i)[2])]:  od:  plot([plt]);

Notice that θ(t) → π as t increases. What motion of the pendulum is represented by this solution?

Exercises: All graphs in the following exercises are to be produced using the plot command, as in the examples above.

1. Derive the equation (1), when there is no drag, i.e., c = 0. [Hint: Recall torque = I·alpha.

2. Display graphs of the solutions sol2 and sol3 on the same coordinate system, for 0 ≤ t ≤ 10. [Hint: Recall the display command in the plots package.]

3. Repeat Exercise 2, but for the initial condition θ(0) = 0, θ′(0) = .5, on the interval 0 ≤ t ≤ 40. Identify which graph is which. [Hint: How long does it take the solution to θ′′ + θ = 0 to complete six periods?]

4. Use conservation of energy to explain why the initial condition θ(0) = 0, θ′(0) = 2 corresponds to the solution of (2) for which the pendulum converges to the vertically upright position, as t → ∞.

Do either Exercise 5 or Exercise 6; extra credit for doing both.

5. Repeat Exercise 2, but add a solution sol4 of the higher-order Taylor approximation to (2),

(3)  \frac{d^{2}θ}{dt^{2}} + θ − \frac{1}{6}θ^{3} + \frac{1}{120}θ^{5} = 0.

Comment on the result. Try also the approximation with the θ^{3} term included but the θ^{5} term omitted; why does this case give trouble?

6. Solve (1) numerically with g/l = 1 and c = .5, using the initial condition θ(0) = 0, θ′(0) = 2. Plot the solution in the θ, θ′ plane. Repeat for different values of θ′(0) until you obtain the solution that corresponds to the pendulum converging to the vertically upright position.