School of Mathematical Sciences

University of Technology, Sydney

35383 High Performance Computing – Autumn 2013

Assignment 1

Due date: Friday 19 April, by noon (12:00 PM)

This assignment is to be completed as an individual task (not as a group work).

1 Purpose of the assignment

The purpose of the assignment is to develop a program (Euler and mid-point methods) for

solving differential equations and apply the methods to the motion of a pendulum. In particular,

you will use your code to answer the following questions:

• If the pendulum is released at an initial angle θ 0 from the vertical, what is the angle at

any time t?

• What is the period of pendulum oscillations, i.e., the time it takes a pendulum to complete

one swing?

Θ

L

m

m g

m g cos Θ

m g sin Θ

Fig. 1: Schematic of the pendulum problem.

2 The pendulum, a nonlinear differential equation

Consider a pendulum with mass m at the end of a rod of length L attached to a fixed pivot

(see Fig. 1). The pendulum normally undergoes a damping because of the friction at the pivot

and the air resistance. But in order to simplify the problem we will assume that the damping is

March 25, 2013

1

negligible (undamped pendulum). The angle between the vertical and the pendulum at time t

is denoted θ(t). The angular equation of motion of the pendulum can be derived from Newton’s

law for circular motion:

−mg L sin(θ(t)) = mL 2

d 2 θ(t)

dt 2

=⇒

d 2 θ(t)

dt 2

+

g

L

sin(θ(t)) = 0. (1)

The pendulum is displaced to an initial angle θ 0 and released at time t = 0 (initial angular

velocity = 0), i.e., we have the initial conditions:

θ(0) = θ 0 (angle of release at t = 0), (2)

dθ(t)

dt

?

t=0

= 0 (angular speed at t = 0). (3)

In Eq. (1), g is the gravitational acceleration constant near the surface of the earth

g = 9.81m/s 2 . (4)

The equation (1) is a nonlinear second order differential equation (with respect to the angle

function θ) and it is not easy to determine its solutions analytically, and so some approximation

techniques are needed in order to solve the problem.

2.1 Approximate solution for small oscillations

For small oscillation angles θ, by using the approximation

sin(θ) ≈ θ, (5)

we can write Eq. (1) as a linear differential equation of second order

d 2 θ(t)

dt 2

+

g

L θ(t)

= 0. (6)

The general solution of the differential equation (6) is

θ(t) = c 1 cos

(

t

√

g

L

)

+ c 2 sin

(

t

√

g

L

)

. (7)

This is a periodic function with period T

T =

2π

√ L/g = 2π

√

L

g

(8)

and the period does not depend on the initial angle θ 0 . Thus for small swing, the period of the

pendulum is almost independent of the maximum angular displacement.

2

3 Numerical methods for the solution of initial value

differential equations

We are interested in the numerical solution of an initial value problem of the form:

dθ(t)

dt

= f(t,θ(t)), (9)

θ(0) = θ 0 . (10)

We assume that this problem has a unique solution over an interval [0,t max ].

3.1 Euler’s method

Euler’s method is one of the simplest numerical methods for solving ordinary differential equa-

tions. In order to approximate the solution, we select N +1 equally spaced values of t over the

interval [0,t max ]:

t i = ih, with h =

t max

N

, (11)

for i = 0,...,N.

In the derivation of Euler’s method, the following approximation is used

dθ(t)

dt

?

t=t i

≈

θ(t i + h) − θ(t i )

h

=

θ(t i+1 ) − θ(t i )

h

. (12)

Substitution into the equation

dθ(t)

dt

= f(t,θ(t)) leads to

θ(t i+1 ) − θ(t i )

h

≈ f(t i ,θ(t i )) =⇒ θ(t i+1 ) ≈ θ(t i ) + hf(t i ,θ(t i )). (13)

The approximation θ(t i+1 ) ≈ θ(t i ) + hf(t i ,θ(t i )) motivates Euler’s method:

θ i+1 = θ i + hf(t i ,θ i ), (14)

where θ i is considered as an estimate to θ(t i ). Note that Euler’s method starts with the initial

value θ 0 .

3.2 Mid-point method

The mid-point method is widely considered as an improvement of Euler’s method and takes

the form:

θ i+1 = θ i + hf

(

t i +

h

2

,θ i +

h

2

f(t i ,θ i )

)

. (15)

3

The following presentation of the mathematical expression (15), as a sequence of compact

formulas, is popular especially for the computer implementations:

k 1 = hf(t i ,θ i ), (16)

k 2 = hf

(

t i +

h

2

,θ i +

1

2

k 1

)

, (17)