Fulltext01 PDF

Modelling and Simulation of Conservative
Dynamical Systems by Computer Algebra

Assisted Lagrangian Mechanics
SA115X Degree Project in Vehicle Engineering, First Level
FREDRIK DEHLIN and JACOB ASK OLSSON
Supervisor: Nicholas Apazidis
Department of Mechanics
School of Engineering Sciences
KTH Royal Institute of Technology
Stockholm, Sweden
May 22, 2017

iii
Abstract
Conservative dynamical systems is modelled with Lagrangian mechanics using

MapleTM with the KTH developed plug-in symbolic package Sophia and sim-
ulated using Matlab® . Two double pendulum configurations and an object
in a Keplerian orbit is studied. Motions and phase portraits are analysed, and
numerical verifications of Kepler’s laws are performed. Properties concerning
chaos is determined partly by examining sensitivity to initial conditions and
it is shown that the 2D pendulum exhibits non-periodic behaviour whilst the
3D pendulum exhibits chaotic behaviour. Kepler’s laws are reproduced under
certain assumptions. Finally, the applicability of Lagrangian mechanics when
applied to conservative dynamical systems is evaluated.
iv
See first, think later, then test.

But always see first. Otherwise,
you will only see what you were
expecting. Most scientists forget
that.
- Douglas Adams
Contents
Contents v
1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.3 Problem Formulations . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3.1 Double Pendulum in 2D . . . . . . . . . . . . . . . . . . . . . 2
1.3.2 Double Pendulum in 3D . . . . . . . . . . . . . . . . . . . . . 3
1.3.3 Planet in Keplerian Orbit . . . . . . . . . . . . . . . . . . . . 4
2 Theory 5
2.1 Newtonian Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1 Newton’s Laws of Motion . . . . . . . . . . . . . . . . . . . . 5
2.1.2 Energy, Work and Conservative Fields . . . . . . . . . . . . . 6
2.1.3 Moment of Inertia . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.4 Coordinate Transforms . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Functionals and Variational Calculus . . . . . . . . . . . . . . . . . . 9
2.3 Lagrangian Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.1 Lagrange’s Equations . . . . . . . . . . . . . . . . . . . . . . 13
2.3.2 Hamilton’s Variational Principle . . . . . . . . . . . . . . . . 13
2.3.3 Lagrange’s Equations from Hamilton’s Principle . . . . . . . 14
2.4 Stability of Dynamical Systems . . . . . . . . . . . . . . . . . . . . . 14
2.5 Orbital Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.5.1 Kepler’s Laws . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.5.2 Mathematics of Ellipses . . . . . . . . . . . . . . . . . . . . . 17
2.5.3 Alternative Derivation of the Radius Vector of an Ellipse . . 22
3 Method 25
3.1 Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 General Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.3 Double Pendulum in 2D . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.1 Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.2 Mechanical Analysis . . . . . . . . . . . . . . . . . . . . . . . 28
v
vi CONTENTS

3.4.1 Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.5 Planet in Keplerian Orbit . . . . . . . . . . . . . . . . . . . . . . . . 33
3.5.1 Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4 Results 37
4.1.1 Trajectory Plots . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.1.2 Phase Portraits . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.2.1 Trajectory Plot . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.2.2 Phase Portraits . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.3.1 Kepler’s First Law . . . . . . . . . . . . . . . . . . . . . . . . 45
4.3.2 Kepler’s Second Law . . . . . . . . . . . . . . . . . . . . . . . 46
4.3.3 Kepler’s Third Law . . . . . . . . . . . . . . . . . . . . . . . . 50
4.3.4 Phase Portraits . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5 Discussion 53
6 Conclusions 57
6.1 Modelling with Lagrangian Mechanics . . . . . . . . . . . . . . . . . 57
6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Bibliography 59
A Maple code 61
A.1 Double Pendulum in 2D . . . . . . . . . . . . . . . . . . . . . 61
A.2 Double Pendulum in 3D . . . . . . . . . . . . . . . . . . . . . 63
A.3 Planet in Keplerian orbit . . . . . . . . . . . . . . . . . . . . 65
B Matlab code 67
B.1 Planet in Keplerian orbit . . . . . . . . . . . . . . . . . . . . 67
B.1.1 twobody.m . . . . . . . . . . . . . . . . . . . . . . . 67
B.1.2 areal_velocity.m . . . . . . . . . . . . . . . . . . . . 72
B.1.3 eccentric_anomaly.m . . . . . . . . . . . . . . . . . 74
Chapter 1
Introduction
1.1 Background
Analytical mechanics is a branch of classical mechanics which deals with new formu-
lations of classical mechanics, where the Lagrangian- and the Hamiltonian formu-
lation are among the most famous ones. This branch of mechanics evolved during
the latter part of the 18th century when mathematicians and physicists discovered
that the Newtonian formulation of classical mechanics did not suffice to solve cer-
tain types of mechanical problems. A new formulation of classical mechanics were
needed and in 1788 one of the most famous ones saw the light of day when the
French-Italian mathematician Joseph-Louis Lagrange published Mécanique analy-
tique containing what is nowadays known as the Lagrangian formulation of classical
mechanics.
With the implementation of Calculus of variations and Hamilton’s principle,
Lagrangian mechanics evolved into a powerful method which can be used to analyse
complex mechanical systems, especially when affected by conservative forces. With
the Lagrangian formulation, an approach based on energy instead of forces paved
the way for revolutionizing theories. Systems containing vast galaxies or microscopic
quantum particles can thus be analysed.
Along with the rise of computers during the middle part of the 20th century,
both previous and newly developed mathematical methods were utilized to take
advantage of the increase in computing power now at hand. The development of
analytical mechanics together with the usage of computers allowed for analysing
virtually any conceivable mechanical system.
1.2 Purpose
The objectives of this thesis are to model and analyse complex mechanical systems
by means of the Lagrangian formulation of classical mechanics with the aid of com-
puter assisted algebra, as well as to perform numerical verification of Kepler’s three
laws of planetary motion.
1
2 CHAPTER 1. INTRODUCTION
1.3 Problem Formulations

This Section introduces the mechanical problems studied symbolically and numeri-
cally in this thesis.
1.3.1 Double Pendulum in 2D

Figure 1.1 shows the rendered Double Pendulum in 2D.
q1
Figure 1.1: Rendered double pendulum in 2D.
The definition of the problem as formulated by Apazidis [1, p. 27]:
’Problem 3.29
Consider a double pendulum consisting of two bars OA and AB of masses m1 and

m2 and lengths `1 and `2 respectively that are free to rotate in the vertical plane
according to the figure. Choose the following numerical values of the parameters
m1 = 3 kg, m2 = 1 kg, `1 = `2 = 1 m and calculate and plot the trajectory
of end B by means of the Sophia and Graphics packages. Choose the following
initial conditions q1 (0) = 1, q̇1 (0) = 0, q2 (0) = 1.5, q̇2 (0) = 0. Calculate then
the trajectory of the point B for slightly different initial conditions and show the
sensitive dependence of the motion of the system on initial conditions by comparing
the two trajectories.’
The sensitivity analysis will be performed by comparing trajectory plots as well as

phase portraits for two cases with different initial conditions.
1.3. PROBLEM FORMULATIONS 3
1.3.2 Double Pendulum in 3D

Figure 1.2 shows the rendered Double Pendulum in 3D.
Figure 1.2: Rendered double pendulum in 3D.
The definition of the problem as formulated by Apazidis [1, p. 26]:
’Problem 3.27
The rigid rod OA of mass m and length ` is pivoted at O so that it may rotate
freely about a horizontal axis through O. A thin rod AB of negligible mass is
mounted at A by means of a ball bearing and may rotate freely about A. A particle
of mass m is attached to the rod at B. Choose some appropriate numerical values
of the parameters and calculate and plot the trajectory of the particle motion for
some initial conditions by means of the Sophia and Graphics packages.’
Apart from presenting the trajectory plot of the particle, various phase portraits
will also be presented and analysed.
4 CHAPTER 1. INTRODUCTION
1.3.3 Planet in Keplerian Orbit

In this problem Kepler’s three laws of planetary motion will be verified by a combina-
tion of numerical and analytical methods. Figure 1.3 shows the two-body planetary
system studied in this problem where the sun is fixed and the planet travels fric-
tionless in the sun’s conservative gravitational field. In Figure 1.3 a, b corresponds
to the semi-major and semi-minor axis respectively whereas c denotes the distance
from the ellipse’s centre to its foci.
b
q1
a c
Figure 1.3: Rendered view of a planet in an elliptical orbit around a star.
Phase portraits for this system will also be analysed and discussed.
Chapter 2
Theory
This chapter introduces the models and underlying theory of classical and analytical
mechanics implemented in this thesis.
2.1 Newtonian Mechanics

To understand why the Lagrangian formulation of classical mechanics is such a
powerful tool, one first must understand the fundamentals of Newtonian mechan-
ics, namely: Newton’s three laws of motion, the relationship between kinetic and
potential energy as well as conservation of energy. In this section, some of these
topics will be covered briefly as an introduction to classical mechanics.
2.1.1 Newton’s Laws of Motion

In 1687 Isaac Newton formulated his three laws of motion which according to
Scheck [2, p. 2] can be translated into something resembling this:
I. ’Every body continues in its state of rest or of uniform rectilinear motion,

except if it is compelled by forces acting on it to change that state.
II. The change of motion is proportional to the applied force and takes place in
the direction of the straight line along which the force acts.
III. To every action there is an equal and contrary reaction; or, the mutual action
of any two bodies are always equal and oppositely directed along the same
straight line.’
These laws are the cornerstones of classical mechanics and by using them the equa-
tions of motion can be formulated for a given non-relativistic macroscopic system.
In general, it would be impossible obtain an analytical solution to the resulting
time-dependent differential equation or system of differential equations since they
are more than likely non-linear.
5
6 CHAPTER 2. THEORY
2.1.2 Energy, Work and Conservative Fields

In classical mechanics, the total mechanical energy of a system can be divided into
two parts, namely kinetic and potential energy. The kinetic energy for a particle is
defined as
m(v · v)
T = , (2.1)
2
whereas the potential energy V (r) is a potential of a conservative function F defined
as Z r
V (r) = − F · dr . (2.2)
r1
To continue the study of energy, the next thing to consider is the relationship
between the two kinds of energies in a situation where only conservative forces acts
on the system. Conservative- forces and forcefields are of great interest in many
disciplines of physics and the underlying theory is found in the mathematical branch
of vector calculus. A conservative vector field is defined as the vector field F where
the line integral around an arbitrarily chosen closed curve C equals to zero [3, p. 28]
I
F · dr = 0 . (2.3)
C
As shown by Matthews [3, Theorem 3.1, p. 51] it is possible to relate the conservative
vector field F to a scalar field V (r) by using the gradient and formulating the
following equation
F = ∇V . (2.4)
Thus by solving (2.4) for V and combining (2.3) with (2.4) the following expression
can be derived Z r2
∇V · dr = V (r2 ) − V (r1 ) , (2.5)
r1
which is of great importance in many physical applications. Equation (2.5) implies
that if you can find the scalar potential V to any vector field you can calculate
the curve integral along an arbitrarily chosen path by evaluating the difference in
potential at the curve’s endpoints.
This concept is applied to physics by using the definition of mechanical work [4,
(9.2), p. 250] Z r2
U1−2 = F · dr (2.6)
r1
and if F is a conservative force field, the work done by that field on a test particle
moving along an arbitrary path C can be calculated with (2.5) as
U1−2 = V (r2 ) − V (r1 ) , (2.7)
where V satisfies (2.4). The relationship between mechanical work and potential
energy is by definition [4, (9.11), p. 253] V1−2 = −U1−2 which for a conservative
force field combined with (2.7) yield the following expression
U1−2 = V1 − V2 . (2.8)
2.1. NEWTONIAN MECHANICS 7
This is a very powerful result since one only have to solve (2.4) to determine if any
force field is conservative. Once a force field is known to be conservative, (2.8) can
be used to determine how much the potential energy has changed from one point to
another in the field. The fact that (2.8) only tells how large the change in potential
is allows you to define your zero point arbitrarily.
To connect the two components of a particles total energy the law of kinetic
energy [4, (9.12), p. 257] is used along with (2.8) to formulate an expression for the
total mechanical energy E as
T +V =E. (2.9)
Equation (2.9) provides a tool to relate the potential energy to the kinetic energy
as long as the particle moves in an conservative force field. Though, this is not the
entire story, when defining the kinetic energy (2.1) the moving object was viewed as
a particle i.e. all of its mass distributed in a single point. In reality when studying
solid objects one might have to consider that the distribution of mass is spread out
in space. This change of physical model results in amongst other things that the
expression for kinetic energy changes, as derived by Apazidis [5, pp. 215-220, 238-
240] from
m(v · v) m(vG · vG ) ω T IG ω
Tparticle = to Ttotal = TG + Trot = + , (2.10)
2 2 2
where vG is the velocity of the object’s centre of mass, ω is the object’s angular
velocity vector and IG is the inertia tensor with respect to the centre of mass G.
2.1.3 Moment of Inertia

In this section, an introduction will be given as to what moment of inertia is and
how it can be computed. As seen in (2.10) the moment of inertia is important to
consider when the rotation of a rigid body is included in the mathematical model,
in contrast to when only a particle is considered. The conclusion from this line of
reasoning is that the moment of inertia has something to do with the distribution
of mass in the rigid body. Since moment of inertia is a tensor quantity it is possible
to divide it into its components. When considering a Cartesian coordinate system
centred in the objects centre of mass those components are as follows [5, p. 156],
Z
Ix = (y 2 + z 2 )dm
ZΩ
Iy = (x2 + z 2 )dm (2.11)
ZΩ
Iz = (x2 + y 2 )dm ,
Ω
where Ω is the rigid body’s volume and dm is an infinitesimal mass element. These
equations are sufficient if the object rotates around one of the coordinate axis, but
in a more generalized case that requirement is not necessarily met.
8 CHAPTER 2. THEORY
In the case of an object rotating around an arbitrary axis Q it can be shown

[5, pp. 214-215] that not only will the components Ix , Iy and Iz in (2.11) appear,
but there will also be a contribution called products of inertia which is defined
accordingly,
Z
Ixy = − xydm
ZΩ
Iyz = − yzdm (2.12)
ZΩ
Izx = − zxdm .
Ω
It can be noted from (2.12) that Ixy = Iyx , Iyz = Izy and Izx = Ixz which together
with (2.11) can be written in a more compact form as the Inertia tensor, IO ,
 
Ixx Ixy Ixz
IO = Iyx Iyy Iyz  . (2.13)
 
Izx Izy Izz
Parallel axis theorem

The parallel axis theorem can be used as a method to transfer the moment of
inertia expressed in a coordinate system centred in the objects centre of mass to
a coordinate system located elsewhere in the object. As to why one would like to
use this theorem can be one of many reasons, but most likely it has to do with the
fact that objects might not rotate around its centre of mass but instead around an
arbitrary axis.
As derived by Apazidis [5, pp. 161-162], it is possible to show that the parallel
axis theorem satisfy the following equation,
IQ = IG + md2 (2.14)
q
where d = x2G + yG 2 is the distance between the new axis of rotation Q and the
centre of mass. There exists an equivalent expression [5, p. 218] for the products of
inertia which is defined as
Ixy = Ix0 y0 − mxG yG , (2.15)
where Ix0 y0 is the product of inertia expressed in the centre of mass and xG , yG
is the distance from the centre of mass to the new coordinate system in x- and y-
direction, respectively.
2.1.4 Coordinate Transforms

An example of where coordinate transforms could be required are when considering
a particle moving in a specific pattern in an arbitrary coordinate system. For a
2.2. FUNCTIONALS AND VARIATIONAL CALCULUS 9
complete understanding of the particles motion it might be necessary to study its

trajectory in a fixed Newtonian system instead of the arbitrary one. To express a
position vector in the fixed system it would be necessary to know the relative motion
as well as the relative inclination between the systems. In the following section the
mathematical theory of this process will be covered briefly.
Rotation of coordinate systems

Consider a Cartesian coordinate system S = (ê1 , ê2 , ê3 ) where an arbitrary vector
v can be written as v = (x1 , x2 , x3 )T . To express v in another Cartesian system
S 0 = (ê01 , ê02 , ê03 ) which is rotated an angle α, β or γ around the ê1 , ê2 or ê3
axis respectively it is necessary to utilize one of the following rotation matrices [6,
pp. 140-141],
1 0 0
 
T1 (α) = 0 cos(α) − sin(α)

 
0 sin(α) cos(α)
cos(β) 0 − sin(β)
 
T2 (β) =  0 1 0 (2.16)
 

sin(β) 0 cos(β)
cos(γ) sin(γ) 0
 
T3 (γ) = − sin(γ) cos(γ) 0 .

 
0 0 1
In the simplest case with only a rotation around one of the coordinate axis e.g.
around ê1 with the angle α, v can be expressed as
v 0 = T1 (α)v (2.17)
in the new coordinate system S 0 .

The transformation matrices (2.16) can also be used to describe an arbitrary
rotation of S 0 by multiplication of corresponding transformation matrices. E.g. if
S 0 is rotated an angle γ around the ê3 -axis and an angle β around the ê2 -axis, v
can be expressed in S 0 as
v 0 = T2 (β)T3 (γ)v . (2.18)
2.2 Functionals and Variational Calculus

A functional is defined as a mapping from a function space to the real or complex
numbers, i.e. a functional takes functions as its input arguments and returns a
scalar. Many, but not all functionals may be expressed in integral form as shown
below. Throughout this text a functional will be denoted with its arguments in
square brackets F [x], whereas parentheses will be used for the argument of a func-
tion.
10 CHAPTER 2. THEORY
A fundamental problem in variational calculus consists of finding a real function

ϕ(x) of a real variable x such that the functional F [ϕ] assumes an extremum. Let
Z x2
d
F [ϕ] = f (ϕ(x), ϕ0 (x), x)dx , ϕ0 (x) ≡ ϕ(x) (2.19)
x1 dx
be this functional with f a given function of ϕ and its derivative, and x1 and x2
are fixed arbitrary endpoints. If ϕ1 = ϕ(x1 ) and ϕ2 = ϕ(x2 ) the problem is to
determine the function ϕ(x) with these endpoints, that will make the functional
F [ϕ] an extremum. Figure 2.1 shows the function ϕ(x) contained within a set of
similar curves.
By defining
Z x2
F () = F [ϕ + η(x)] = f (ϕ(x) + η(x), ϕ0 (x) + η 0 (x), x)dx (2.20)
x1
as a small variation in the argument of the functional, where η and ϕ is considered

fixed and η(x1 ) = 0 = η(x2 ). The function F () is a function of and its extremum
can be determined by differentiating with respect to and setting the derivative
equal to zero. Differentiating results in
Z x2
∂f ∂f 0

0
F () = η(x) + η (x) dx (2.21)
x1 ∂ϕ ∂ϕ0
and by using partial integration on the second term
Z x2
∂f d ∂f

0
F () = − η(x)dx , (2.22)
x1 ∂ϕ dx ∂ϕ0
where the previously stated conditions on η(x) at the endpoints has been used.
φ2
φ(x)
φ1
x1 x2 x
Figure 2.1: The curve ϕ(x) makes the functional F [ϕ] an extremum and is contained
within a set of curves with the same endpoints.
2.3. LAGRANGIAN MECHANICS 11
The expression within brackets in (2.22) is called the functional derivative of f with
respect to ϕ and in order for ϕ to be a stationary function this expression must be
equal to zero for all variations η(x). This results in the differential equation
∂f d ∂f

− =0, (2.23)
∂ϕ dx ∂ϕ0
which is called Euler-Lagrange’s differential equation of variational calculus in one
dimension.
2.3 Lagrangian Mechanics

Lagrangian mechanics is a fundamental part of general mechanics and allows for the
study of more complex and abstract problems than those that can be constructed
and solved with Newtonian mechanics. Lagrangian mechanics introduces gener-
alized coordinates that are connected to the constraints of the system as well as
concepts such as configuration space.
Generalized Coordinates
A set of coordinates that are independent and include the spatial constraints of the
system can be labeled as generalized coordinates. This concept can be exemplified
by a particle moving on the surface of a sphere with radius R, as in Figure 2.2.
y
z x
Figure 2.2: The xy-plane of a sphere showing the particle moving in a circular
motion on a circle with radius R from the origin.
The constraints of this system can be expressed as x2 + y 2 + z 2 = R2 , in the depen-

dent variables {x, y, z}. By introducing spherical coordinates {r, θ ϕ} the motion
of the particle can be represented by q1 = θ and q2 = ϕ which are independent
variables and generalized coordinates.
As stated by Scheck [2, p. 89]:
’In general, the set of 3n space coordinates of the n-particle system will be replaced
by a set of (3n − Λ) generalized coordinates, as such
r1 , r2 , ..., rn → q1 , q2 , ..., qf , f = 1, 2, ..., 3n − Λ ’ , (2.24)
where r denotes the space coordinates, q the generalized coordinates and f the
degrees of freedom as a function of the number of particles n and the number of
constraints Λ.
In the example illustrated in Figure 2.2, n = 1 since one particle is considered
and Λ = 1 because the particle is constrained to the surface of the sphere. This
consequently results in the two generalized coordinates, q1 = θ and q2 = ϕ.
Configuration Space
By specifying the configuration of a system by the use of generalized coordinates

(i.e. specifying the position of the particles), the configuration space represents all
possible configurations of the system. The dimension of the configuration space
is equal to the minimum number of coordinates needed to completely specify a
configuration. The dimension can also be interpreted as the number of degrees of
freedom (compare with f in (2.24)), if the constraints of the system are holonomic 1 .
A coordinate function χ can be used to parameterize the configuration space and
the configuration path γ is a mapping from time to configuration-space points [7].
Coordinate Path
The coordinate path is essentially a mapping from time to tuples2 of generalized

coordinates. For each instant of time t, the coordinate path can be represented by
q1 (t)
 
 q2 (t) 
 
q(t) =  .   , (2.25)
 .. 

qn (t)
which describes the motion of the system. The coordinate path corresponds to the
configuration path by q(t) = χ(γ(t)) where χ is the coordinate function and γ is
the configuration path.[7]
1
A holonomic constraint is such that it may be expressed in the following form:
f (q1 , q2 , ..., qn , t) = 0, where {q1 , ..., qn } are the generalized coordinates of the system. [2]
2
A tuple is a finite ordered list of elements.
2.3. LAGRANGIAN MECHANICS 13
2.3.1 Lagrange’s Equations

By examining the virtual displacements δri of a system of n-particles one can obtain
an equation expressing d’Alembert’s principle of virtual displacements:
n
(Ki − ṗi ) · δri = 0 , (2.26)
X
i=1
where Ki is the real, dynamic force and ṗi is the momentum of the i-th particle.
For a holonomically constrained system, the virtual displacements can be ex-
pressed as a function of the generalized coordinates. By defining the virtual displace-
ments in terms of generalized coordinates, the quantity δqk will be independent, in
contrast to δri from (2.26) which are dependent. As shown by Scheck [2, pp. 92-94],
replacing δri in (2.26) by δqk the following set of equations is obtained:
d ∂T ∂T

− = Qk , k = 1, 2, ..., 3n − Λ , (2.27)
dt ∂ q˙k ∂qk
where Qk is a quantity called generalized forces, T is the kinetic energy of the sys-
tem and Λ is the number of constraints.
Furthermore Sheck shows [2, p. 94] that if the generalized forces Qk are conser-
vative, they can be expressed as
∂
Qk = −V , (2.28)
∂qk
where V is the potential energy of the system, a function of the generalized coordi-
nates and time t. By introducing the Lagrangian function, L as
L(qk , q˙k , t) = T (qk , q˙k ) − V (qk , t) , (2.29)

then (2.27) takes the form
d ∂L ∂L

− =0, (2.30)
dt ∂ q˙k ∂qk
which are called Lagrange’s equations [2].
2.3.2 Hamilton’s Variational Principle

If it is shown that the system evolves along a path that satisfies (2.29), it is possible
to define the action integral in the same form as (2.19) in Section 2.2, namely
Z t2
S[q] = L(q(t), q̇(t), t)dt , (2.31)
t1
where L is the Lagrangian of the system. This restatement of the principle of
stationary action in terms of energies is called Hamilton’s variational principle. By
defining q(t) as a solution to the equations of motion that assumes the boundary
values q(t1 ) = x1 and q(t2 ) = x2 , the curve q(t) will be such that the action integral
(2.31) will assume an extremum [7].
2.3.3 Lagrange’s Equations from Hamilton’s Principle

As stated in Section 2.3.2, the integral curve q(t) must be defined such that the
action integral (2.31) assumes an extremum. This integral curve must satisfy the
Lagrange’s equations:
∂L d ∂L

− =0, k = 1, 2, ..., 3n − Λ . (2.32)
∂qk dt ∂ q˙k
This statement is proven analogously to (2.23) in Section 2.2 [2].
2.4 Stability of Dynamical Systems

In an effort to thoroughly describe the behaviour of a dynamical system, the con-
cept of stability and chaos is a useful matter of discussion. There exists a wide
variety of analytical tools in this subject. This chapter is an introduction to one of
these, namely the Phase Portrait and how it can represent the sensitivity to initial
conditions of the system.
Phase Portrait
A phase portrait is a geometric representation of the trajectories of a dynamical
system when viewed in the phase plane.
Consider a single pendulum as in Figure 2.3.
L
g
=q
m
Figure 2.3: Single pendulum of a mass m and length L, with the angle θ as a
generalized coordinate q and g the acceleration due to gravity.
In the absence of friction and no external driving force, the pendulum is in a con-
servative force field due to g being the only acting acceleration. Since the system is
conservative, the total mechanical energy E defined as in (2.9) will be constant.
By defining the displacement angle θ from the y-axis as the generalized coordi-
nate q, it can be shown that the phase portrait for this system looks schematically
as in Figure 2.4 [8, pp. 168-170].
2.4. STABILITY OF DYNAMICAL SYSTEMS 15
.
q
Figure 2.4: Phase portrait of the single pendulum for different energy levels.
In Figure 2.4 the trajectories have been presented for different energy levels as con-
tours. The generalized coordinate q and its derivative with respect to time q̇ is
represented on the x-axis and y-axis respectively.
The physical interpretation of the phase portrait can be summarized as follows
[8, p. 170].
• The centre point (0, 0) is the stable equilibrium point where the pendulum is
at rest hanging vertically.
• Small orbits around the centre represent small oscillations about equilibrium,
i.e. the pendulum moves in a back-and-forth manner exhibiting periodic mo-
tion.
• As E increases, the pendulum will reach a state where the mass moves in a
continuous circular motion around the centre.
It is worth noting that even though the phase portrait generally has a strong con-
nection to the energy of the system, it is not necessary to generate the portrait as
a contour-plot for different energy levels. A set of data for q(t) and q̇(t) where t is
the time, is sufficient to generate the phase portraits and draw conclusions about
the stability of the system.
An attractor is defined as a region in the state space (i.e. a set of numerical
values) toward which the system tends to evolve. In contrast, a repellor is a region
that tends to repel the trajectories. If a system does not seem to behave in any
periodic manner and the number of degrees of freedom is greater than two, it can
usually be labelled as chaotic.
2.5 Orbital Mechanics

Orbital mechanics is the principal study of the motion of bodies moving around
other bodies under the influence of gravity.
2.5.1 Kepler’s Laws

Johannes Kepler determined a mathematical description of planetary motion and
formulated three laws.
As stated by Hintz3 [9, p. 5]:
’Kepler’s Laws are:
I. The orbit of each planet is an ellipse with the sun at a focus.
II. The line joining the planet to the sun sweeps out equal areas inside the ellipse
in equal time intervals. Therefore, the velocity at closest approach is greater
than the velocity at the furthest distance from the sun.
III. The square of the period of a planet is proportional to the cube of its mean
distance from the sun.’
Kepler’s first and second laws are illustrated in Figure 2.5.
Figure 2.5: Keplerian motion of a planet around a star.
The Second Law indicates that vp > va which can be realized by stating that in
order for the two grey-colored area segments in Figure 2.5 to be equal, the velocity
of the orbiting planet near the perihelion must be greater than the velocity near the
aphelion.
3
The prevalent definition of Kepler’s Third Law is normally formulated as ’The square of the
period of a planet is proportional to the cube of the semi-major axis of its orbit.’
2.5. ORBITAL MECHANICS 17
2.5.2 Mathematics of Ellipses

There are some important concepts as well as mathematical theories concerning
ellipses which needs to be covered before any more rigorous analysis of Kepler’s
laws can be performed. Consider the ellipse in Figure 2.6 with a polar coordinate
system centred in its right focus. The ellipse is characterised by a semi-major and
semi-minor axis of length a and b respectively. In many mathematical texts it is
possible to find expressions for the radius as well as the sector area, though they are
usually derived with respect to a coordinate system located in the ellipses’ centre.
In orbital mechanics and especially when proving Kepler’s laws, it is necessary to
describe the radius as well as the sector area with respect to a focus centred coor-
dinate system.
Figure 2.6: Ellipse with a polar coordinate system in its right focus.
To derive an expression for the radius r(θ) and the sector area A(θ), one would
first have to consider an ellipse with a Cartesian coordinate system centred in the
ellipses’ centre. The eccentricity e is defined [10] as
s
a

e≡ 1− , (2.33)
b
from which it is possible to derive the following equations relating the semi-major
and semi-minor axis, as well as the distance c between the foci and the centre to
the eccentricity,
p
b = a 1 − e2 (2.34)
c = ae . (2.35)
Inserting a polar coordinate system in the right focus as shown in Figure 2.6 results
in the coordinate transform defined as
x = r(θ) cos(θ) + c
(
(2.36)
y = r(θ) sin(θ) .
An ellipse centred in origo is defined as
2 2
x y
+ =1 (2.37)
a b
and by inserting (2.36) into (2.37) an expression for the ellipse with respect to the
polar coordinate system can be formulated as
2 2
r cos(θ) + c r sin(θ)

+ = 1. (2.38)
a b
Expanding (2.38) and inserting relations (2.34) and (2.35), the following expression
is obtained
a2 (1 − e2 )(r2 cos2 (θ) + 2aer cos(θ) + a2 e2 ) + a2 r2 (1 − cos2 (θ)) = a4 (1 − e2 ) , (2.39)
which simplifies into

r = ±(er cos(θ) − a(1 − e2 )) . (2.40)
Setting e = 0 it is possible to determine which sign in (2.40) is the correct one.
When e = 0, r = ±a. Since a > 0 one would have to use the negative sign in
equation (2.40) seeing that the positive sign would give a negative radius. By using
that result it is possible to solve (2.40) for the radius r(θ), resulting in
a(1 − e2 )
r(θ) = . (2.41)
1 + e cos(θ)
To calculate the area A(θ), the area integral I is evaluated over the entire region
Ω covered by the ellipse,
ZZ ZZ
I= dA = {dA = dxdy} = dxdy , (2.42)
Ω Ω
and by changing into polar coordinates
r(θ)2
Z θ Z r(θ) Z θ
A(θ) = I = rdrdθ = dθ . (2.43)
0 0 0 2
If (2.41) and (2.43) are combined it is possible to calculate the radius as well as the
area from the right focus for an arbitrary ellipse with eccentricity e. The case when
the polar coordinate system is placed in the right focus, one would have to redo the
calculations from equation (2.36) down to equation (2.41). Though in that case it
is important to remember that the x component is expressed as x = r(θ) cos(θ) − c

instead of x = r(θ) cos(θ) + c which is used in this thesis.
For many applications in celestial mechanics as well as in astrodynamics it is of
great importance to be able to calculate the time to a specific point in the elliptical
orbit. These calculations can be performed using Kepler’s equation [9, p. 49],
s
a3
tP = (ε − e sin(ε)) , (2.44)
µ
where tP is the time elapsed from the periapsis to position P in the orbit, a is
the semi-major axis, e is the eccentricity and ε is the eccentric anomaly. The
gravitational parameter µ is defined [9, p. 26] as
µ = G(m1 + m2 ) , (2.45)
where G is the universal constant of gravitation and m1 , m2 is the mass of the two
particles.
The derivation of (2.44) will not be shown here. For an in depth derivation,
Chapter 2 in Orbital Mechanics and Astrodynamics by Hintz, G. [9, pp. 23-53] is
highly recommended.
As can be seen in (2.44) the eccentric anomaly ε is required to determine the
time to a specific point in the orbit. When the true anomaly θ is known it is possible
to derive an expression for ε as a function of θ. The true- and eccentric anomalies
is defined in Figure 2.7 below.
Figure 2.7: Defining the eccentric- and true anomaly.
Consider the ellipse in Figure 2.7 with a polar coordinate system in its right
focus and a Cartesian coordinate system in its centre. The vector F P 4 which runs
from focus F to the arbitrary point P on the ellipse is thus written as
!
r(θ) cos(θ) + ae
FP = , (2.46)
r(θ) sin(θ)
where r(θ) is the radius defined in (2.41), a is the semi-major axis and e is the
eccentricity defined in (2.33).
Next step is to multiply the y-component by a/b as a way to stretch the ellipse
in the y-direction, thus making it into a circle. Introduce the vector CQ which is
written as
4
A vector will in this section be defined in bold where the letters define the point it begin and
end at respectively.
!
r(θ) cos(θ) + ae
CQ = = {Insert (2.41)} =
b r(θ) sin(θ)
a
!
a (1 − e2 )√
cos(θ) + e + e2 cos(θ)
= . (2.47)
1 + e cos(θ) 1 − e2 sin(θ)
Using the result in (2.47) it is possible to formulate an expression for the eccentric
anomaly as
√
sin(ε) 1 − e2 sin(θ)
tan(ε) = = . (2.48)
cos(ε) (1 − e2 ) cos(θ) + e + e2 cos(θ)
By introduction of the Weierstrass substitution which according to Spivak is
’The worlds sneakiest substitution’ [11, pp. 382-383], the following relations hold
true.
Let 
2t
 sin(x) =
1 + t2

x

t = tan , then 2 . (2.49)
2 cos(x) = 1 − t


1 + t2
Equation (2.48) is evaluated with the Weierstrass substitution (2.49), beginning
with the left hand side of the equal sign and moving on to the right hand side after
that. This gives
2t
sin(ε) 2 tan 2ε

1+t2
= = , (2.50)
cos(ε) 1−t2 1 − tan2 2ε
1+t2
√ √
1 − e2 sin(θ) 1 − e2 1+t
2t
2
= 2 , (2.51)
(1 − e2 ) cos(θ) + e + e2 cos(θ) 1−t2
(1 − e ) 1+t2 + e + e2 1−t
2
2
1+t
which after some algebra becomes
q
2 1−e
1+e tan θ
2
... = . (2.52)
1− 1−e
1+e tan2 θ
2
Inserting (2.50) and (2.52) into (2.48)

q
2 tan 2ε 2 1−e
tan θ

1+e 2
ε =
, (2.53)
1 − tan 2
2

1− 1−e
1+e tan2 θ
2
and after identification, an expression for the eccentricity can be formulated as

s
ε 1−e θ

2 tan =2 tan , (2.54)
2 1+e 2
which after simplification becomes

s
1−e
!
θ
ε(θ) = arctan tan . (2.55)
1+e 2
By inserting (2.55) into Kepler’s equation (2.44) an expression for the time it takes
from periapsis to a specific point in the orbit can be calculated as a function of the
true anomaly, s
a3
tp (θ) = (ε(θ) − e sin(ε(θ))) . (2.56)
µ
2.5.3 Alternative Derivation of the Radius Vector of an Ellipse

The aforementioned derivation of equation (2.41) might be unintuitive for a reader
not familiar with the Weierstrass substitution and the eccentric anomaly. Therefore,
an alternative derivation using only properties of the ellipse is presented below.
Consider the ellipse in Figure 2.8 with a semi-major axis a and a distance from
its centre to a focus c.
y
r2
r1
focus 1 θ
f
x
focus 2
Figure 2.8: Defining the eccentric- and true anomaly.
The vector f in Figure 2.8 is defined as
f ≡ 2c êx , (2.57)
which enables the following equation to be formulated
r2 = f + r1 . (2.58)
Using the following trait of an ellipse
r1 + r2 = 2a (2.59)
and taking the square of the magnitude of (2.58)
||r2 ||2 = ||f + r1 ||2 ⇐⇒ r22 = 4c2 + 4cr1 cos(θ) + r12 , (2.60)
which after utilization of the conjugate rule and (2.59), (2.61) can be simplified into
2 (a − r1 ) 2a = 4c (c + r1 cos(θ)) . (2.61)
By reformulating (2.35) as
c
e= (2.62)
a
and after inserting (2.62) into a rewritten form of (2.61) it results in
r1 c c r1 r1 r1

1− = + cos(θ) ⇐⇒ 1 − = e2 + e cos(θ) . (2.63)
a a a a a a
After some further simplifications (2.63) can be written as
a 1 − e2

r1 (θ) = , (2.64)
1 + e cos(θ)
which is exactly the same expression as in (2.41).

Chapter 3
Method
This chapter describes the methods used to model and solve the three mechanical
problems presented in Section 1.3.
3.1 Software
To model, solve and present the mechanical problems various software were used.
This section briefly describes these programs and their respective area of application.
Maple and Sophia-package

Maple provides the user with a symbolical and numerical computing environment.
The theory presented in Chapter 2 that deals with modelling the mechanical system
is implemented in Maple.
Due to the extensive number of equations needed in the computation of a me-
chanical system, a package called ’Sophia’ is loaded and used in the Maple envi-
ronment. Sophia is a symbolic package that facilitates formulation of equations of
motion for a mechanical system in symbolic form. In this thesis it is done by means
of Lagrangian formulation. Sophia allows for a simplified method of constructing
the problem and consequently a more compressed code. The various commands
exclusive to the Sophia-package that has been utilized are summarized and briefly
described below.
- dependsTime declares a variable to be time dependent
- &rot establishes the relation between two coordinate systems
- &ev evaluates vector in a specific coordinate system
- toTimeFunction substitution set used to differentiate with respect to time
- &dt differential tool with respect to time
By using Maple and Sophia the rigid body problems studied in this thesis were
constructed and solved efficiently for the generalized coordinates at different time
25
26 CHAPTER 3. METHOD
steps. By exporting these values to an Excel-sheet, further analysis with different

software were made possible.
The Maple code can be seen in Appendix A.
Matlab
Matlab is primarily a program intended for numerical computing and visualization
through graphs and plots. By importing the solutions from Maple into Matlab
for the different mechanical problems, various graphs were created to present data.
Using Matlab’s animation functions, video sequences of the analysed mechanical
systems were recorded.
Blender
Blender is a free and open source 3D animation software developed and maintained
by the Blender Foundation. Blender is primarily used in this thesis as a tool to
visualize trajectories and motions of the studied problems, but it is also used to
render some of the images in this report.
3.2 General Procedure

To analyse the mechanical problems an algorithmic approach with regard to mod-
elling and solving were utilized. This method can be summarized as such:
I. Define the generalized coordinates, constants and frame relations by examining

the geometry of the mechanical system.
II. Construct the Lagrangian L = T − V , where T is the kinetic energy of the

system and V the potential energy.
III. Formulate the Lagrange equations of motion.
IV. Solve the generalized coordinates and their respective derivatives numerically
using a Fehlberg fourth-fifth order Runge-Kutta method [12].
When the abovementioned steps have been completed, the problem can be further
analysed and visualized using Matlab and Blender.
The following sections in this chapter presents the method applied to each prob-
lem in a more thorough manner.
3.3. DOUBLE PENDULUM IN 2D 27
3.3 Double Pendulum in 2D

This section describes the geometry and mechanics of the Double Pendulum when
analysed in two dimensions.
Figure 3.1 shows the rendered double pendulum.
O
z
x
q1

Figure 3.1: Rendered double pendulum in 2D with inserted coordinate systems N ,

A and B.
3.3.1 Modelling
The double pendulum consists of two bars OA and AB of masses m1 and m2 and
lengths `1 and `2 respectively. The bars can rotate freely without friction in the
vertical xy-plane. The system exists in a conservative force field with standard
gravitational acceleration g.
A fixed reference frame N as well as two rotating Cartesian coordinate systems
A and B are defined as in Figure 3.1. Both A and B rotates around their respective
third-axis, i.e. around êA B
3 and ê3 . The generalized coordinates are q1 and q2 ,
where q1 is defined as the angle between bar OA and the y-axis of N . The second
coordinate q2 are defined as the angle between bar AB and the êA
2 -axis.
3.3.2 Mechanical Analysis

In accordance with the proposed procedure in Section 3.2 the second step consists
of determining the kinetic and potential energy of the system in order to define the
Lagrangian L.
Equation (2.29) states that L is defined as the difference between the kinetic
and the potential energy of the system, and by applying the theory explained in
Section 2.1.2, the Lagrangian is calculated as follows.
Kinetic Energy
The system consists of two rigid objects, namely bar OA and bar AB. In order to
calculate the kinetic energy for the bars as proposed by (2.10), the velocity vector
of the centre of masses are required.
The velocity vector for bar OA is
h i
vA
1 =
1
2 `1 q̇1 0 0 (3.1)
and for bar AB
h i
vB
2 = 2 `2 (q̇1
1
+ q̇2 ) + `1 cos (q2 )q̇1 `1 sin (q2 )q̇1 0 , (3.2)
where (3.1) is defined in the A-frame and (3.2) in the B-frame. This is denoted by
the use of superscripted A and B.
The moments of inertia for the bars are
1
I1 = m1 `21 (3.3)
3
and
1
I2 =
m2 `22 , (3.4)
12
for bar OA and AB respectively. Equation (3.3) has been derived using the parallel
axis theorem (2.14) and it represents the moment of inertia about the end of the
bar located in O. Equation (3.4) is the moment of inertia about the centre of mass
for bar AB, which is located in the middle of the bar. Note that both bars are
presumed to have negligible thickness and uniformly distributed mass.
Since the endpoint of the bar is fixed in O, the contribution to the kinetic energy
of the system will only depend on the rotational motion and I1 is determined as
in (3.3). The bar AB will however experience both translational and rotational
motion. By using (2.10) the kinetic energy for the bars is
1
T1 = I1 q̇12 (3.5)
2
and
1 1
T2 = I2 (q̇1 + q̇2 )2 + m2 (v B B
2 · v2 ) . (3.6)
2 2
Potential Energy
The reference level for the conservative force field is set to be the x-axis of the N -
frame and the potential energy for the two bars are therefore defined by the vertical
displacement of the centre of masses with regard to this level.
By analysing the geometry of the system, the potential energy for bar OA is
calculated to
1
V1 = − m1 `1 g cos (q1 ) (3.7)
2
and for bar AB to
1

V2 = −m2 g `1 cos (q1 ) + `2 cos (q1 ) cos (q2 ) . (3.8)
2
Lagrangian
In order to determine the Lagrangian L for this system, (2.29) is utilized by sum-
mation of the energy contributions for the two bars.
The Lagrangian for the double pendulum in 2D is
L = T1 + T2 − (V1 + V2 ) , (3.9)
where T1 , T2 , V1 and V2 are defined by (3.5)-(3.8).
Parameters, Lagrange equations and solving

After assigning values to the parameters m1 , m2 , `1 , `2 and g as well as initial
conditions for the generalized coordinates, the Lagrange equations can be defined for
the mechanical system. By applying a numerical solver, the generalized coordinates
and their derivatives are obtained at discrete times.
For this specific problem, a set of two equations are derived which can be found
in Appendix A.1.

This section describes the geometry and mechanics of the Double Pendulum when
analysed in three dimensions.
Figure 3.2 shows the rendered double pendulum.
Figure 3.2: Rendered double pendulum in 3D with inserted coordinate systems N ,

A, B and C.
3.4.1 Modelling
The double pendulum consists of two bars, bar OA with mass m and bar AB with
negligible mass. At the end of bar AB at point B, a particle with mass m is located.
The length of both bars are `. The bar OA may rotate freely without friction about
a horizontal axis through O and bar AB may rotate freely about A. The system
exists in a conservative force field with standard gravitational acceleration g.
A fixed reference frame N as well as three rotating frames A, B and C are
Cartesian coordinate systems and are defined as in Figure 3.2. The origin of N is
located in point O, A is located at a distance of `/2 between point O and A and B

is located at point A. The fourth coordinate system C, is located at point B inside
the particle.
The generalized coordinates are q1 , q2 and q3 , where q1 is defined as the angle
between bar OA and the z-axis of N , and q2 is the angle between the êA 3 - and
êB
3 -axis. The third generalized coordinate q 3 is defined as the angle between the êB
3-
C
and ê3 -axis.

In accordance with the proposed procedure in Section 3.2 the second step consists
of determining the kinetic and potential energy of the system in order to define the
Lagrangian L.
Equation (2.29) states that L is defined as the difference between the kinetic
and the potential energy of the system, and by applying the theory explained in
Section 2.1.2, the Lagrangian is calculated as follows.
Kinetic Energy
The system consists of two rigid objects namely the bars OA and AB. In order to
calculate and present the trajectory for the particle located in B, the motion of the
two bars must be considered.
The bar OA will only experience rotational motion about the fixed point O.
The kinetic energy for bar OA can therefore be calculated as
1
T1 = I1 q̇12 , (3.10)
2
by applying (2.10). The moment of inertia-term I1 in (3.10), has been derived using
the parallel axis theorem (2.14) resulting in
1
I1 = m`2 . (3.11)
3
Since the bar AB is considered massless, the second contribution to the kinetic
energy of the system consists of the motion of the particle in B, which naturally
is dependent on the motion of the two bars. As calculated in Appendix A.2, the
velocity vector of the particle is
h i
v C2 = `(cos (q2 )q̇1 + cos (q3 )q̇1 + q̇3 ) `(sin (q2 ) sin (q3 )q̇1 + cos (q3 )q̇2 ) −` sin (q3 )q̇1
(3.12)
which is a vector defined in frame C.
By applying (2.10), the kinetic energy T2 for the particle can be derived as
1
T2 = m(v C2 · v C2 ) . (3.13)
2
Potential Energy
The reference level for the conservative force field is set to be the y-axis of the N -
frame and the potential energy for the bar OA and the particle in B are therefore
defined by the vertical displacement of the centre of masses with regard to this level.
The potential energy for the bar OA is
1
V1 = − m`g cos (q1 ) , (3.14)
2
which can be realized by inspecting the geometry in Figure 3.2.
In order to define the potential energy for the particle, the position vector to the
point B 1 with regard to the fixed N -frame is required. Using the Sophia-package,
this position vector is derived to
−` [cos (q2 ) sin (q1 ) cos (q3 ) + cos (q1 ) sin (q3 ) + sin (q1 )]
 
rN C = ` sin (q2 ) cos (q3 )  , (3.15)

 
−` [cos (q2 ) cos (q1 ) cos (q3 ) − sin (q1 ) sin (q3 ) + cos (q1 )]
which is the displacement of the particle with regards to the reference level. The
potential energy for the particle is therefore

V2 = mg r N C · êN
z . (3.16)
Lagrangian
In order to determine the Lagrangian L for this system, (2.29) is utilized by sum-
mation of the energy contributions for the bar OA and the particle.
The Lagrangian for the double pendulum in 3D is
L = T1 + T2 − (V1 + V2 ) , (3.17)
where T1 , T2 , V1 and V2 are defined by (3.10), (3.13), (3.14) and (3.16) respectively.

After assigning values to the parameters m, ` and g as well as initial conditions for
the generalized coordinates, the Lagrange equations can be defined for the mechan-
ical system. By applying a numerical solver, the generalized coordinates and their
derivatives are obtained at discrete times.
For this specific problem, a set of three equations are derived which can be found
in Appendix A.2.
1
Note that the point B coincides with the origin of the C-frame.
3.5. PLANET IN KEPLERIAN ORBIT 33
3.5 Planet in Keplerian Orbit

This section describes the geometry and mechanics of the orbital motion of the
planet. Figure 3.3 below shows a rendered view of the planetary system.
Figure 3.3: Rendered view of a planet in an elliptical orbit around a star.
3.5.1 Modelling
When studying planetary motion it is possible to formulate different kinds of math-
ematical models depending on what the purpose of the study is. In this case an
analysis of a planet orbiting around a significant heavier object, e.g. a star, is per-
formed. Since both objects exerts an equal but opposite gravitational force on each
other (in accordance with Newtons third law motion) and the mass of a star [13] is
many orders of magnitude greater than the mass of a planet, approximating the star
as stationary in space is reasonable. The system is also considered to only consist
of particles since the effect of the rigid body’s rotation plays a minor part in the
objects motion on a planetary scale.
To allow for easier visualization as well as interpretation of the result from the
forthcoming analysis the systems’ parameters have been chosen to other values than
the astronomical ones. Mass of the star and mass of the orbiting planet is denoted
m1 and m2 respectively. The particle corresponding to the planet orbits frictionless
throughout space and the sole force acting on the two particles is gravity described
by Newton’s law of universal gravitation [2, p. 8],
mi mk
Fki = −G r̂ik , (3.18)
kri − rk k2
where Fki is the gravitational force between the two particles (i and k) and G is the
universal gravitational constant. The distance between the two particles is denoted
kri − rk k whereas r̂ik is the unit vector from particle i to k and mi , mk is the
particles respective mass.
A fixed reference polar coordinate system has been inserted with the star fixated
in its centre as shown in Figure 3.3. The generalized coordinates are q1 and q2 , where
q1 is defined as the radius to the planet and q2 is defined as the angle θ measured
from the line between the star and periapsis.

Kinetic Energy
The system consists of two particles where one of them, namely the star, is kept
stationary. As discussed in Section 3.5.1 the objects are considered to be particles,
resulting in that equation (2.1) is used to calculate the kinetic energy of the system.
The velocity vector for the second particle located in point A when expressed in
the A-frame is derived to h i
v1A = q̇1 q1 q̇2 0 , (3.19)
resulting in the expression for the kinetic energy written as
1 1
T = m2 (v1A · v1A ) = m2 (q12 q̇22 + q̇12 ) . (3.20)
2 2
Potential Energy
The potential energy function for a particle in a gravitational field is derived by
using the relations between mechanical work and potential energy from Section
2.1.2 and by solving (2.4) for the scalar field V with (3.18) as prevailing vector field.
An expression for the gravitational potential function is derived as
mi mk
V (r) = −G , (3.21)
r
where r = kri − rk k is the distance between particles i and k. By using (3.21)
with the aid of the generalized coordinates an expression for the total amount of
potential energy in the two body system can be formulated as
mi mk
V = −G . (3.22)
q1
Lagrangian
The Lagrangian of this system is formulated by insertion of (3.20) and (3.22) into
(2.29) resulting in
L=T −V . (3.23)

After the Lagrangian has been formulated and numerical values for the constants
G, m1 , m2 as well as for the initial has been assigned, the Lagrange equations can
be formulated. By using the Runge-Kutta-Fehlberg numerical differential equation
solver of order four and five, values for the generalized coordinates and their respec-
tive time derivatives is obtained for discrete points in time. The actual Lagrange
equations derived when solving this problem can be found in Appendix A.3.
Chapter 4
Results
This chapter presents the results of the analysed problems described in Section 1.3.

Table 4.1 shows the parameters used to define the properties of the system.
Parameters m1 = 3 kg, m2 = 1 kg, l1 = 1 m, l2 = 1 m, g = 9.8 m/s2

Table 4.1: Parameters for the double pendulum in 2D.
Table 4.2 shows the initial conditions for the generalized coordinates and their
derivatives for the two different cases A and B.
Initial Conditions
Case q1 (0) [rad] q̇1 (0) [rad/s] q2 (0) [rad] q̇2 (0) [rad/s]
A 1 0 1.5 0
B 1.1 0 1.4 0
Table 4.2: Initial conditions for the 2D pendulum in case A and case B.
Note that the parameters in Table 4.1 are defining properties in both cases.
The different initial conditions will alter the motion of the double pendulum as
time progresses. The resulting differences between the cases are presented later in
this section in the form of trajectory plots and phase portraits.
37
38 CHAPTER 4. RESULTS
4.1.1 Trajectory Plots

Figure 4.1 shows the trajectory of the point B for case A in 4.1a and case B in
4.1b.
Trajectory Plot for point B Trajectory Plot for point B
0.5 0.5
0 0
-0.5 -0.5
-1 -1
-1.5 -1.5
-2 -2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2
(a) Trajectory for point B, case A. (b) Trajectory for point B, case B.
Figure 4.1: Trajectory plots for point B for t : 0 → 10 s, for the different cases A
and B.
The trajectories are generated for a time interval of ten seconds.
4.1.2 Phase Portraits

Figure 4.2 shows the phase portrait for the generalized coordinate q1 for case A in
4.2a and for case B in 4.2b. The portraits are generated for a time interval of 100
seconds.
Phase Portrait: q 1(t) Phase Portrait: q 1(t)

6 6
4 4
2 2
0 0
-2 -2
-4 -4
-6 -6
-1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5
(a) Generalized coordinate q1 , case A. (b) Generalized coordinate q1 , case B.
Figure 4.2: Phase portraits for the generalized coordinate q1 (t) for t : 0 → 100 s,
for the different cases A and B.
The phase portraits for the generalized coordinate q2 are presented in Figure 4.3,
where 4.3a and 4.3b describes case A and B respectively.

20 20
15 15
10 10
5 5
0 0
-5 -5
-10 -10
-15 -15
-20 -20
-70 -60 -50 -40 -30 -20 -10 0 10 20 30 -70 -60 -50 -40 -30 -20 -10 0 10 20 30
(a) Generalized coordinate q2 , case A. (b) Generalized coordinate q2 , case B.
Figure 4.3: Phase portraits for the generalized coordinate q2 (t) for t : 0 → 100 s,
for the different cases A and B.
Note that the x-axis (generalized coordinate as a function of time) in the phase
portraits is not restricted to an interval of ±2π. The values of around −65 rad in
Figure 4.3b can be interpreted as a rotational motion of the bar AB in a clockwise
manner for a large number of consecutive revolutions.
Figure 4.4 shows a sequence in which both cases A (blue) and B (green) are pre-
sented simultaneously.
Geometric Con guration Phase Portrait: q1(t) Phase Portrait: q2(t)

2 5
15
Angular velocity [rad/s]

Spatial coordinate [m]
4
1.5
3 10
1
2
5
0.5 1
(a) 0
0
-1
0
-0.5 -5
-2
-1 -3
-10
-1.5 -4
-5 -15
-2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -1.5 -1 -0.5 0 0.5 1 1.5 -2 0 2 4 6 8 10 12 14
2 5
15

4
1.5
3 10
1
2
5
0.5 1
(b) 0
0
-1
0
-0.5 -5
-2
-1 -3
-10
-1.5 -4
-5 -15
-2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -1.5 -1 -0.5 0 0.5 1 1.5 -2 0 2 4 6 8 10 12 14
2 5
15
4
1.5
3 10
1
2
5
0.5 1
(c) 0
0
-1
0
-0.5 -5
-2
-1 -3
-10
-1.5 -4
-5 -15
-2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -1.5 -1 -0.5 0 0.5 1 1.5 -2 0 2 4 6 8 10 12 14
2 5
15
4
1.5
3 10
1
2
5
0.5 1
(d) 0
0
-1
0
-0.5 -5
-2
-1 -3
-10
-1.5 -4
-5 -15
-2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -1.5 -1 -0.5 0 0.5 1 1.5 -2 0 2 4 6 8 10 12 14
2 5
15

4
1.5
3 10
1
2
5
0.5 1
(e) 0
0
-1
0
-0.5 -5
-2
-1 -3
-10
-1.5 -4
-5 -15
-2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -1.5 -1 -0.5 0 0.5 1 1.5 -2 0 2 4 6 8 10 12 14
Spatial coordinate [m] Angle [rad] Angle [rad]
Figure 4.4: Geometric configuration, Phase portraits for q1 (t) and Phase portraits
for q2 (t), in left, middle and right column respectively. Time after start of simula-
tion: (a) ∆t = 4.2 s; (b) ∆t = 4.5 s; (c) ∆t = 4.9 s; (d) ∆t = 5.1 s; (e) ∆t = 6.5
s.
The left column of subfigures shows the simulated geometric configuration of both
double pendulums, and the middle and right column shows the phase portraits
for q1 and q2 respectively. The rows indicate various times after the start of the
simulation sequence.
Figure 4.5 shows the phase portraits for both coordinates and cases, A in blue
and B in green.
Phase Portrait: q1(t) Phase Portrait: q2(t)

5
15
4
3 10
2
5
1
0 0
-1
-5
-2
-3
-10
-4
-5 -15
-1.5 -1 -0.5 0 0.5 1 1.5 -2 0 2 4 6 8 10 12 14
(a) Generalized coordinate q1 . (b) Generalized coordinate q2 .
Figure 4.5: Phase portraits for the generalized coordinate q1 (t) and q2 (t)
for t : 0 → 10 s, for both cases A and B.

In this section the results obtained when analysing the double pendulum in 3D is
presented.
Initial conditions as well as the parameters defining the system can be found in
Table 4.3 below.
Parameters m = 1 kg, l = 1 m, g = 9.8 m/s2

q1 (0) = 1 rad, q̇1 (0) = 1 rad/s,
Initial Conditions q2 (0) = 1 rad, q̇2 (0) = 1 rad/s,
q3 (0) = 1 rad, q̇3 (0) = 2 rad/s
Table 4.3: Parameters and initial conditions for the double pendulum in 3D.
4.2.1 Trajectory Plot

When analysing the motion of the double pendulum in 3D it is of interest to plot
the trajectory of the particle and utilizing it as a tool to further increase the under-
standing of the particle’s motion. In Figure 4.6 below the trajectory is plotted for
a time interval of ten seconds.
Trajectory Plot for Particle
0.5
-0.5
-1
-1.5
-2
1
0.5 2
0 1
0
-0.5
-1
-1 -2
Figure 4.6: Trajectory of the particle with initial conditions and parameters defined
in Table 4.3.
From Figure 4.6 it can be noted that with the aforementioned initial conditions and
parameters the particle will follow an irregular and virtually chaotic trajectory.

Phase portraits for the three generalized coordinates q1 , q2 and q3 is generated for
a time interval of 10 seconds and is presented in Figure 4.7 below.

15 300
10 200
5 100
0 0
-5 -100
-10 -200
-15 -300
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 0 5 10 15 20 25 30 35
Phase Portrait: q 3(t)

25
20
15
10
-5
-10
-15
-20
-25
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
(c) Generalized coordinate q3 .
Figure 4.7: Phase portraits for the generalized coordinates for t : 0 → 10 s.

Phase portraits for the three generalized coordinates q1 , q2 and q3 is generated for
a time interval of 100 seconds and is presented in Figure 4.8 below.

15 1600
1400
10
1200
1000
5
800
0 600
400
-5
200
0
-10
-200
-15 -400
-3 -2 -1 0 1 2 3 -100 -80 -60 -40 -20 0 20 40
Phase Portrait: q3(t)

30
20
10
-10
-20
-30
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
(c) Generalized coordinate q3 .
Figure 4.8: Phase portraits for the generalized coordinates for t : 0 → 100 s.
Note that the angle which correspond to the generalized coordinate q2 in Figure 4.8b
is not restricted to an interval of ±2π. Instead an angle which is greater or lesser
than ±2π can be interpreted as consecutive revolutions in the same direction. To
exemplify, q2 reaches values around −100 rad which corresponds to a large number
of consecutive clockwise revolutions.

In this section the results from the analysis of a planet in a Keplerian orbit will
be presented along with numerical proofs of the validity of Kepler’s three laws of
planetary motion. The parameters and initial conditions used can be found in Table
4.4 below.
Parameters m1 = 250 kg, m2 = 1 kg, G = 0.5 m3 /(kg · s2 )
Initial Conditions q1 (0) = 0.34 m, q̇1 (0) = 0 m/s,
q2 (0) = 0 rad, q̇2 (0) = 75 rad/s
Table 4.4: Parameters and initial conditions for the planet in a Keplerian orbit.
4.3.1 Kepler’s First Law

Kepler’s first law states that
The orbit of each planet is an ellipse with the sun at a focus.
To verify this, different initial conditions have been used in the numerical equation
solver described in Section 3.5.2 producing different orbits for the planet. The
equation for an ellipse with its centre in (h, k) is given by
2 2
x−h y−k

+ =1 (4.1)
a b
where a is the semi-major axis and b is the semi-minor axis. By inserting the coor-
dinates of an arbitrarily selected point in the orbit into (4.1) along with respective
semi-major and semi-minor axis lengths it should equal to one in case the orbit is
an ellipse.
Results obtained when entering different initial conditions in the numerical solver
can be found in Table 4.5 below.
2 2
Initial conditions (a, b) (h, k) x−h
a + y−k
b
q1 (0) = 0.34, q2 (0) = 0 a = 1.4698 h = −1.1298 min = 1, max = 1
q̇1 (0) = 0, q̇2 (0) = 75 b = 0.9401 k=0 mean = 1, median =1
q1 (0) = 0.40, q2 (0) = 0 a = 2.5510 h = −2.1510 min = 1, max = 1
q̇1 (0) = 0, q̇2 (0) = 60 b = 1.3714 k=0 mean = 1, median =1
q1 (0) = 0.50, q2 (0) = 0 a = 0.4545 h = 0.0455 min = 1, max = 1
q̇1 (0) = 0, q̇2 (0) = 30 b = 0.4523 k=0 mean = 1, median =1
q1 (0) = 0.20, q2 (0) = 0 a = 2.3656 h = −2.1656 min = 1, max = 1
q̇1 (0) = 0, q̇2 (0) = 173 b = 0.9520 k=0 mean = 1, median =1
Table 4.5: Verification of Kepler’s first law for different initial conditions with pa-
rameter values m1 = 250, m2 = 1 and G = 0.5. Analysis is performed on all of the
calculated x- and y- coordinates (≈ 10000 points).
As can be seen in the rightmost column in Table 4.5, when testing different initial
conditions the planet follows an elliptical orbit since the left hand side of (4.1)
always equals to one. Kepler’s first law can therefore be considered to be valid.
4.3.2 Kepler’s Second Law

Kepler’s second law states that
The line joining the planet to the sun sweeps out equal areas
inside the ellipse in equal time intervals. Therefore, the velocity
at closest approach is greater than the velocity at the furthest
distance from the sun.
Three different methods was used to verify the validity of Kepler’s second law. Two
of the three methods involved the use of Kepler’s equation (2.44) and the sector
area of an ellipse (2.43) derived in Section 2.5.2 whilst the third method used the
sector velocity.
I. Numerical Approach with MATLAB

The main idea of this method was to numerically evaluate the sector area for each
true anomaly (θ) obtained from the numerical solution of the Lagrange equations.
This first method can be summarized as follows.
By inserting (2.41) into (2.43) the swept area can be numerically evaluated as a
function of θ using matlab’s integral command. It is thus possible to calculate
the swept area for each discrete time step. This process is repeated iteratively until
the swept area equals a predetermined value at which the corresponding elapsed
time is recorded. The same process is repeated starting from the aphelion and the
elapsed time until reaching the predetermined area is recorded.
Figure 4.9 shows that the predetermined area was selected to π/12 area units.
t 2 = 0.0608 s t 1 = 0.0605 s
- Orbiting Planet - Sun
Figure 4.9: Numerical verification of the validity of Kepler’s second law of planetary
motion.
As seen in Figure 4.9 the radius vector swept A1 during ∆t1 = 0.0605 s, and A2
during ∆t2 = 0.0608 s.
II. Kepler’s Equation and Sector Area

The second method used to verify Kepler’s second law of planetary motion was
to calculate both the sector area for each true anomaly with (2.43) as well as the
elapsed time from the periapsis with Kepler’s equation (2.56) for each true anomaly.
If Kepler’s second law holds true the result when plotting swept area as a function
of elapsed time should be a first order equation with constant derivative.
When performing the aforementioned method, area swept by the radius vector as a
function of elapsed time since periapsis is obtained and presented in Figure 4.10.
Swept area as a function of elapsed time

5
4.5
3.5
Swept area [area unit]
2.5
1.5
0.5 Calculated data

Reference line: A = 4.34t + 0.3
0
0 0.2 0.4 0.6 0.8 1
Elapsed time since periapsis [s]
Figure 4.10: Area swept by the radius vector as a function of elapsed time since
periapsis.
It is of interest to note a few important observations from this figure. Firstly, as

expected from Kepler’s second law the swept area obeys a linear equation. Secondly,
the derivative is constant with a value of ~4.34 area units per second.
Once the swept area is calculated as a function of elapsed time it is possible to
verify Kepler’s second law to an even greater extent. By linearising the derivative
of the swept area with respect to time, the sector velocity can be written as
dAi ∆Ai Ai+1 − Ai

≈ = for i = 1, . . . , n . (4.2)
dt ∆ti ti+1 − ti
Evaluation of (4.2) for each discrete time results in the sector velocity as a function
of elapsed time since periapsis, as shown in Figure 4.11 below.
Sector velocity as a function of elapsed time since periapsis

4.344
4.3439
Sector velocity [(area unit)/s]
4.3438
4.3437
4.3436
4.3435
4.3434
4.3433
0 0.2 0.4 0.6 0.8 1
Figure 4.11: Sector velocity of the radius vector as a function of elapsed time since
periapsis.
Figure 4.11 confirms the conclusions drawn from Figure 4.10, i.e. Kepler’s second
law of planetary motion is valid.
III. Sector Velocity

The third method used to verify Kepler’s second law of planetary motion was the
use of an expression for the sector velocity derived by Hintz [9, p. 52], defined as
dA r2
= θ̇ , (4.3)
dt 2
which by using the generalized coordinates defined in Section 3.5.1 can be written
as
dA q2
= 1 q̇2 . (4.4)
dt 2
By inserting the numerically computed values for the generalized coordinates ob-
tained from the Lagrange equations into (4.4), the sector velocity can be calculated
as a function of the discrete time steps.
When calculating the sector velocity using (4.4) Figure 4.12 below is obtained.
Sector velocity as a function of elapsed time since periapsis

4.34
4.339
4.338
Sector velocity [(area unit)/s]
4.337
4.336
4.335
4.334
4.333
4.332
4.331
4.33
0 0.2 0.4 0.6 0.8 1
Figure 4.12: Sector velocity of the radius vector as a function of elapsed time since
periapsis.
This method results in a sector velocity of ~4.335 area units per second. The sector
velocity has now been utilized to verify Kepler’s second law of planetary motion.
4.3.3 Kepler’s Third Law

Kepler’s third law states that
The square of the period of a planet is proportional to the cube

of the semi-major axis of its orbit.
To verify the validity of Kepler’s third law of planetary motion one first consider
Kepler’s equation (2.44). By inserting an eccentric anomaly of E(θ) = 2π, (2.44)
simplifies to s
a3
tp (θ = 2π) = 2π , (4.5)
µ
which can be rewritten as

2
tp

µ = a3 ⇒ τ 2 ∝ a3 (4.6)
2π
where τ is the period i.e. time for one orbit and a is the semi-major axis. Through
this argumentation it is confirmed that Kepler’s third law of planetary motion is
valid.

Phase portraits for the two generalized coordinates q1 and q2 is generated for a time
interval of 2.5 seconds and is presented in Figure 4.13 below.

15 80
70
10
60
5
50
0 40
30
-5
20
-10
10
-15 0
0 0.5 1 1.5 2 2.5 3 0 2 4 6 8 10 12 14 16
Figure 4.13: Phase portraits for the generalized coordinates for t : 0 → 2.5 s.
It can be noted from Figure 4.13a that the generalized coordinate q1 which corre-
sponds to the radius as well as q2 which corresponds to the true anomaly exhibits
periodic patterns. This is expected from a planet moving in a stable Keplerian
orbit.
Chapter 5
Discussion
This chapter contains discussions concerning the results presented in Chapter 4.

The trajectory plots for point B in Figure 4.1 shows distinct differences in motion
between the two cases A and B after only a few seconds. This can be interpreted
as an indication that the system is highly sensitive to initial conditions.
The phase portraits in Figure 4.2 for the first generalized coordinate q1 (the
angle between the first bar OA and the vertical y-axis of the N system) are not
periodic since no closed-loop orbits are seemingly repeated. The phase space is
becoming increasingly covered as time progresses, which is a further indication of
a non-periodic behaviour. It is noticeable that there are differences between the
two phase spaces corresponding to case A and to case B in Figure 4.2a and 4.2b
respectively. This is emphasized more clearly in Figure 4.5a where the phase portrait
was generated during a time of ten seconds.
Figure 4.3 shows the phase portraits for the second generalized coordinate q2
for the two cases A and B. The portraits suggest a quasi-periodic behaviour of the
coordinate. The most prominent difference between the two cases is the location of
the attractors in Figure 4.3a compared to those in Figure 4.3b. In case A there seem
to be more subspaces in the phase space that exhibits an attractive nature, compared
to case B. The locations of the attractors are also different. The attractors for case
B is located between 0 to 20 rad, whereas the ones apparent in case A are located
between −30 and −10 rad as well as around 12 rad. It is also noticeable that q2
reaches larger negative values in case B, to around −65 rad, which means that the
bar AB completes a greater amount of consecutive clockwise revolutions in case B
in comparison to case A.
The sequence shown in Figure 4.4 highlights some of the discrepancies between
the two cases during the start of the simulation. At t = 0, the point A of the
pendulum in case B is located slightly above the same point in case A, whilst the
point B is located roughly in the same location for both cases. This results in a
53
54 CHAPTER 5. DISCUSSION
higher energy level for case B, which can also be noted by inspecting the phase
portrait for q1 at row (a) in Figure 4.4. The case B (green) trajectory follows the
case A (blue) rather closely, aside from being located farther away from the origin
and thus generating a larger phase space which can be seen both in Figure 4.2 or
Figure 4.5a. By further examining the sequence it is possible to pinpoint when the
two cases start to behave differently. From t = 0 to t = 4.2 seconds, i.e. row (a)
in Figure 4.4, both phase portraits trace similar trajectories. Between 4.5 and 4.9
seconds the motion of the two pendulums differ largely, as can be seen in the phase
portraits at rows (b) and (c). The last row, (e), shows the state of the pendulums
after 6.5 seconds and the difference between the cases are distinct. The trajectory
in case A (blue) for q2 seem to be attracted to a region of around 13 rad whilst the
trajectory in case B (green) seem to be in a quasi-periodic orbit, visible in Figure
4.5b.
The sequence in Figure 4.4 indicates that the double pendulum is sensitive to
initial conditions and the effects can be seen after a few seconds.

The motion of the double pendulum in 3D is characterized by abrupt transitions,
i.e. a kind of jerkiness can be observed. In Figure 4.6, this can be seen as sudden
changes in direction of the trajectory. To fully understand this behaviour, one must
examine the actual modelling of the dynamical system as well as the initial condi-
tions used in the simulation.
The bar AB is assumed to have negligible mass and the two generalized coor-
dinates, q2 and q3 is therefore governed by the motion of the particle in point B
with mass m. The initial conditions specified in Table 4.3 results in a large amount
of kinetic energy being introduced into the system. Since the effect of friction is
excluded in this model, the energy of the system will be conserved. These factors
might explain the spikes which can be seen in the phase portrait for q2 in Figure
4.7b and 4.8b, and consequently the jerkiness of the motion.
The phase portraits in Figure 4.7 and 4.8 are not exhibiting any periodic be-
haviour and the phase space for q1 and q3 are becoming increasingly covered as time
progresses.

First and foremost, what needs to be considered when analysing the results obtained
from the study of Kepler’s three laws of planetary motion is that the parameters
used (which can be found in Table 4.4) differs from the astronomical ones. As an
example, the Sun’s actual mass is in the order of magnitude of 1031 kg whereas a
mass of 250 kg have been used in this thesis. The universal gravitational constant
is ten orders of magnitude greater than in reality to compensate for the Sun having
a vastly lower mass in this model.
The main reason as to why these changes in values were made was to minimize
the effects of numerical errors. The floating point precision of Matlab is 16 dec-
imals which means that to remain within acceptable accuracy one would have to
change the Sun’s mass roughly 1014 kg. To circumvent these types of computational
problems occurring when working with values on an astronomical scale, the possi-
bility to change into- and perform the calculations in the astronomical system of
units [14] exists. In this system of units, the base for length is the astronomical unit
(au) which corresponds to the distance from the Sun to Earth, instead of meter as
in the SI-system. Another example is that instead of kg, the IAU-system use one
solar mass as its base for mass.
Figure 4.13 shows the phase portraits for the generalized coordinates q1 and
q2 corresponding to the radius and the true anomaly respectively. Both phase
portraits demonstrate a periodic behaviour of the coordinates which have to be the
case when a planet orbits in a Keplerian orbit. An object which follows e.g. a
hyperbolic trajectory will not exhibit any periodic behaviour since its orbit is not
bound to the larger body. It will thus not obey Kepler’s laws of planetary motion.
Kepler’s first law of planetary motion
The results obtained when verifying Kepler’s first law of planetary motion were
conclusive. The planet followed an elliptical orbit regardless of initial conditions. It
is worth to note that the reason as to why q2 (0) and q̇1 (0) have been kept equal to
zero for all different cases (as shown in Table 4.5) solely has to do with programming
issues. It turned out to be problematic to determine the length of the semi-major
and semi-minor axis from a discrete set of data points when the ellipse was tilted.
Instead it was decided that only initial conditions which kept the major and minor
axis parallel to the coordinate system were to be used.
Kepler’s second law of planetary motion
Kepler’s second law of planetary motion was verified by several different methods,
each with its own strengths and weaknesses. The first method involved an iterative
process where the elapsed area during small changes in time were summed up to a
predetermined value. This process was repeated at a later stage in the orbit and
the two elapsed times were compared. The discrepancy between the two calculated
times as seen in Figure 4.9, are most likely due to numerical errors. Since the input
to (2.43) are discrete values it is possible that the calculated areas ends up either
slightly bigger or smaller than the predetermined value of π/12. This results in a
difference in the calculated time since A1 and A2 may differ in size. To minimize
this error a smaller time step can be used when solving the Lagrange equations.
As shown in Section 4.3.2 the second method used was to plot the swept area as
a function of elapsed time by means of Kepler’s equation. As expected, this process
resulted in a linear function (see Figure 4.10) which after differentiation proved that
56 CHAPTER 5. DISCUSSION
the sector velocity is constant for all times (see Figure 4.11). This method resulted
in a sector velocity of approximately 4.344 area units per second.
The third and final method used to verify Kepler’s second law involved utilizing
an analytically derived expression for the sector velocity. Insertion of the generalized
coordinates computed by Maple into (4.4) resulted in the plot shown in Figure 4.12
with the constant sector velocity of approximately 4.335 area units per second.
A discrepancy in the result between the last two methods can be noted. The
reason for this has to do with the accuracy of the two different methods used. It is
expected that the last method (i.e. the analytical) is the most accurate of the two
since it rests firmly on an analytical foundation whereas the second method requires
more numerical calculations.
Weaknesses of Kepler’s laws

Even though it is remarkable that planetary orbits can be calculated using physical-
and mathematical models originating with Kepler and Newton in the 17th to 18th
century, they are merely approximations. A requirement for Kepler’s laws is that
the object with largest mass remains stationary in one of the foci. It is a quite
good approximation when considering objects where the difference in mass is many
orders of magnitude, but the approximations falters when comparing e.g. Jupiter
to our Sun. Due to Jupiter’s enormous mass [15] only being around three orders of
magnitude less than the Sun’s mass, their point of rotation (i.e. their barycentre)
is located a few thousand kilometres outside of the Sun’s surface [16, p. 199]. The
Sun will consequently orbit a point in space with a radius of approximately 700 000
km thus rendering the approximation of it being stationary invalid.
Another aspect which needs to be considered when using Kepler’s laws is that
they assume an isolated system where the planet is only affected by a single source
of gravity, namely the Sun in this case. In reality this is far from true. To correctly
calculate the trajectories of planets or other objects one has to consider the pertur-
bations caused by all of the other planets’ gravitational pulls. For instance, when
Earth and Jupiter are aligned on opposite sides of the Sun, Jupiter’s gravitational
pull enhances that of the Sun’s whilst when aligned on the same side Jupiter will
cancel out some parts of the Sun’s gravitational pull. This causes an alteration in
Earth’s orbit which is not accounted for by Kepler’s laws.
Finally, the model used in this thesis has its foundation in classical Newtonian
mechanics. This means that relativistic effects are not accounted for. As a con-
sequence, especially when studying objects orbiting close to another object with
enormous mass, the Newtonian formulation of classical mechanics has to be aban-
doned in favour of general relativity.
Chapter 6
Conclusions
This chapter draws conclusions based on the results presented in Chapter 4 and the
discussions presented in Chapter 5.
6.1 Modelling with Lagrangian Mechanics

To model dynamical systems with the use of the Lagrangian formulation of classical
mechanics is a highly effective and direct method. Some might find the underlying
theory of Lagrangian mechanics difficult to comprehend but the concept of mod-
elling by energy is an intuitive approach when analysing more complex mechanical
systems, i.e. systems containing several degrees of freedom.
In order to derive the necessary Lagrange equations that defines a system, one
might discover that the size of the algebraic expressions tend to become substantial
as well as complicated. This is especially true when the system is described by
numerous generalized coordinates. By using sophisticated software (such as Maple
with the Sophia-package, as used in this thesis) the task to algebraically derive these
equations becomes significantly simpler.
Lagrangian mechanics provides a reliable framework, with which a large range of
problems may be analysed. Apart from solving the equations of motion, Lagrangian
mechanics also reveals conserved quantities and their symmetries in a direct manner.
Together with Hamilton’s Variational Principle, Lagrangian mechanics is proven
useful in a variety of branches within theoretical physics, particularly in quantum
mechanics. By modelling and simulating a seemingly simple conservative dynamical
system such as the double pendulum, the complexity of the mechanical motion is
revealed.
57
58 CHAPTER 6. CONCLUSIONS
The Pendulums
A crucial difference between the two pendulums is the additional degree of freedom
introduced in the 3D problem. This has a large influence on, not only the actual
motion of the pendulum but also on the characteristic signature of the system as
can be seen in the phase portraits. While 2D problems with two degrees of freedom
is at most non-periodic, 3D problems with three degrees of freedom may be chaotic.
The fact that a dynamical system can be chaotic and thus sensitive to initial con-
ditions, provides an explanation as to why the complex motion arises. By analysing
the phase portraits, certain behavioural aspects of the system can be interpreted.
In order to compare and qualitatively discuss the differences between the pen-
dulums, additional tools not used in this thesis is required.
Keplerian Orbit
The use of Lagrangian mechanics when analysing Kepler’s three laws of planetary
motion was not necessary in the two-body problem since there exists an analytical
solution. The need for Lagrangian mechanics arises when considering an unre-
stricted three-body problem, though it is also here one of the largest drawbacks of
the Lagrangian formulation arises. In order to model a system with three bodies,
one would have to define three generalized coordinates for each body resulting in a
total of 9 Lagrange equations to solve.
6.2 Future Work

The study of chaos and stability contains a wide variety of concepts including
analysing underlying patterns, feedback loops, fractals, symmetries etc. Theories
concerning attractors and quasi-periodicity can further enhance the understanding
of obtained results. As an addition to phase portraits, the use of Poincaré maps
might provide useful information.
Noethers Theorem can be applied at an early stage of the analysis to state con-
servation laws governing the system.
To increase the accuracy of celestial models, the Hamiltonian formulation of
classical mechanics is advantageous over the Lagrangian formulation since the pos-
sibility to include multiple bodies is more direct.
Bibliography
[1] Apazidis N. Computer Algebra Assisted Mechanics; 2017.
[2] Scheck F. Mechanics : From Newton’s Laws to Deterministic Chaos. Berlin,

Heidelberg : Springer Berlin Heidelberg; 2010.
[3] Matthews PC. Vector Calculus. 1st ed. Springer Undergraduate Mathematics
Series. Springer-Verlag London; 1998. 7th printing 2005.
[4] Apazidis N. Mekanik. 1, Statik och partikeldynamik. 2nd ed. Lund: Lund :
Studentlitteratur; 2013.
[5] Apazidis N. Mekanik. 2, Partikelsystem, stel kropp och analytisk mekanik. 1st
ed. Lund: Lund : Studentlitteratur; 2012.
[6] Arfken GB, Weber HJ, Harris FE. Mathematical Methods for Physicists. 7th ed.
Boston: Academic Press; 2013. Available from: http://www.sciencedirect.
com/science/article/pii/B978012384654900030X.
[7] Sussman GJ, Wisdom J. Structure and Interpretation of Classical Mechanics.

2nd ed. Cambridge, Massachusetts; London, England: MIT Press; 2015.
[8] Strogatz SH. Nonlinear dynamics and chaos : with applications to physics,
biology, chemistry and engineering. 2nd ed. Cambridge, Massachusetts.: West-
View Press; 2000.
[9] Hintz GR. Orbital Mechanics and Astrodynamics - Techniques and Tools for
Space Missions. 1st ed. Springer International Publishing; 2015.
[10] Weisstein EW. "Ellipse.". MathWorld–A Wolfram Web Resource.; 2017.

[Accessed 11th April 2017]. [Online]. Available from: http://mathworld.
wolfram.com/Ellipse.html.
[11] Spivak M. Calculus. 3rd ed. Houston, Texas: Pubish or Perish, Inc.; 1994.
[12] "dsolve/numeric/rkf45". Maplesoft Website;. [Accessed 10th April 2017].

[Online]. Available from: http://www.maplesoft.com/support/help/Maple/
view.aspx?path=dsolve/rkf45.
59
60 BIBLIOGRAPHY
[13] Williams DR. "Sun Fact Sheet". NASA Goddard Space Flight Center; 2016.
[Accessed 19th April 2017]. [Online]. Available from: https://nssdc.gsfc.
nasa.gov/planetary/factsheet/sunfact.html.
[14] IAU. "Measuring the Universe: The IAU and astronomical units". International
Astronomical Union; 2017. [Accessed 9th May 2017]. [Online]. Available from:
https://www.iau.org/public/themes/measuring/.
[15] Williams DR. "Jupiter Fact Sheet". NASA Goddard Space Flight Center; 2016.
[Accessed 11th May 2017]. [Online]. Available from: https://nssdc.gsfc.
nasa.gov/planetary/factsheet/jupiterfact.html.
[16] MacDougall DW. Newton’s Gravity: An Introductory Guide to the Mechanics

of the Universe. 1st ed. Springer-Verlag New York; 2012.
Appendix A
Maple code
In this Appendix the Maple code used can be found.
A.1 Double Pendulum in 2D

> restart;
> with(plots): with(plottools): with(ArrayTools): with(ExcelTools):
> read ‘C:/path/SophiaV6N.txt‘:
> read ‘C:/path/Graphics.txt‘:
> # ––––––– I. Modelling the system ––––––– #
> dependsTime(q1,u1,q2,u2):
> &rot[N,A,3,q1]:
> &rot[A,B,3,q2]:
> r__O := N &ev [0,0,0]:
> r__A := r__O &++ (A &ev [0,-l1/2,0]):
> r__B := r__O &++ (A &ev [0,-l1,0]) &++ (B &ev [0,-l2/2,0]):
> r__P := r__O &++ r__B &++ (B &ev [0,-l2/2,0]):
> # Defining parameter values & initial conditions
> param := {m1=3,m2=1,l1=1,l2=1,g=9.8}:
> InitCond := {q1(0)=1.1,u1(0)=0,q2(0)=1.4,u2(0)=0}:
> # ––––––– II. Mechanical Analysis ––––––– #
> v__1 := simplify(N &fdt r__A):
> v__2 := simplify(N &fdt r__B):
> I__1 := 1/3*m1*l1^2:
> I__2 := 1/12*m2*l2^2:
> T__1 := 1/2*I__1*q1t^2:
> T__2 := 1/2*I__2*(q1t+q2t)^2 + 1/2*m2*(v__2 &o v__2):
> T := simplify(T__1 + T__2):
> V := -1/2*l1*m1*g*cos(q1)-m2*g*(l1*cos(q1)+1/2*l2*cos(q2)*cos(q1)):
> L := T-V:
61
> # ––––––– III. Euler-Lagrange equations ––––––– #
> L := subs(q1t=u1,q2t=u2,L):
> L__q1 := diff(L,q1):
> L__q2 := diff(L,q2):
> L__u1 := diff(L,u1):
> L__u2 := diff(L,u2):
> L__u1_t := &dt L__u1:
> L__u2_t := &dt L__u2:
> eq1 := L__u1_t-L__q1=0:
> eq1 := subs(q1t=u1,q2t=u2,eq1):
> eq2 := L__u2_t-L__q2=0:
> eq2 := subs(q1t=u1,q2t=u2,eq2):
> kde := {q1t=u1,q2t=u2}:
> eq1 := {u1t=solve(eq1,u1t)}:
> eqs := eq1 union eq2 union kde:
> # ––––––– IV. Solver ––––––– #
> eqst := subs(toTimeFunction,eqs):
> eqst := subs(param,eqst):
> ff := dsolve(eqst union InitCond,{q1(t),q2(t),u1(t),u2(t)},type=numeric,
maxfun=500000):
> numberOfPoints := 720:
> integrationTime := 30:
> _plot := odeplot (ff,[[t,q1(t)],[t,q2(t)],[t,u1(t)],[t,u2(t)]],0..
integrationTime , numpoints = numberOfPoints):
> gg := dsolve(eqst union InitCond,{q1(t),q2(t),u1(t),u2(t)},type=numeric,
maxfun=500000,output=Array([seq((1/100*i),i=0..10000)])):
> gg := gg(2,1):
> # Exporting results

> Export(gg,"Problem_3.29_B_tot.xls","Problem_3.29")
A.2 Double Pendulum in 3D
> restart;
> read ‘C:/path/SophiaV6N.txt‘;
> read ‘C:/path/Graphics.txt‘;
> dependsTime(q1,u1,q2,u2,q3,u3):
> &rot[N,A,2,q1]:
> &rot[A,B,1,q2]:
> &rot[B,C,2,q3]:
> r__O := N &ev[0,0,0]:
> r__A := r__O &++ (A &ev[0,0,-l/2]):
> r__B := r__O &++ (A &ev[0,0,-l]) &++ (B &ev[0,0,0]):
> r__C := r__O &++ r__B &++ (C &ev[0,0,-l]):
> r__CN := N&to r__C:
> r__CN := r__CN[1]:
> InitCond := {q1(0)=1,u1(0)=1,q2(0)=1,u2(0)=1,q3(0)=1,u3(0)=2}:
> param := {m=1,l=1,g=9.8}:
> v__A := simplify(N &fdt r__A):
> v__B := simplify(N &fdt r__B):
> v__C := simplify(N &fdt r__C):
> I__1 := 1/3*m*l^2:
> T__1 := 1/2*I__1*q1t^2:
> T__2 := 1/2*m*(v__C &o v__C):
> T := simplify(T__1+T__2):
> V__1 := -l/2*m*g*cos(q1):
> V__2 := m*g*(r__CN[3]):
> V := V__1+V__2:
> L := simplify(T-V):
> L := subs(q1t=u1,q2t=u2,q3t=u3,L):
> L__q1 := diff(L,q1):
> L__q2 := diff(L,q2):
> L__q3 := diff(L,q3):
> L__u1 := diff(L,u1):
> L__u2 := diff(L,u2):
> L__u3 := diff(L,u3):
> L__u1_t := &dt L__u1:
> L__u2_t := &dt L__u2:
> L__u3_t := &dt L__u3:
> eq1 := L__u1_t-L__q1=0:
> eq1 := subs(q1t=u1,q2t=u2,q3t=u3,eq1):
> eq2 := L__u2_t-L__q2=0:
> eq3 := L__u3_t-L__q3=0:
> kde := {q1t=u1,q2t=u2,q3t=u3}:
> eqs := eq1 union eq2 union eq3 union kde:
> # ––––––– IV. Solver ––––––– #
> eqst := subs(toTimeFunction,eqs):
> eqst := subs(param,eqst):
> ff := dsolve(eqst union InitCond,{q1(t),q2(t),q3(t),u1(t),u2(t),u3(t)},
type=numeric,maxfun=500000):
> plotpath_3d(ff,r__C,param,0,10,800,black,0,N):
> gg := dsolve(eqst union InitCond,{q1(t),q2(t),q3(t),u1(t),u2(t),u3(t)},
type=numeric,maxfun=500000,output=Array([seq((1/240*i),i=0..3600)])):
> gg := gg(2,1):

> Export(gg,"Problem_3.27.xls","Problem_3.27"):
A.3 Planet in Keplerian orbit
> restart;
> read ‘C:/path/SophiaV6N.txt‘;
> read ‘C:/path/Graphics.txt‘;
> dependsTime(q1,u1,q2,u2):
> &rot[N,A,3,q2]:
> r__O:= N &ev[0,0,0]:
> r__A:= r__O &++ (A &ev[q1,0,0]):
> v__A:= simplify(N &fdt r__A):
> InitCond:={q1(0)=0.34,u1(0)=0,q2(0)=0,u2(0)=75}:
> param:={m = 1, M = 250, G = 5*10^(-1)}:
> T := 1/2*m*(v__A &o v__A):
> V:= -m*G*M/q1:
> L:=T-V:
> L:=subs(q1t=u1,q2t=u2,L):
> L__q1:=diff(L,q1):
> L__q2:=diff(L,q2):
> L__u1:=diff(L,u1):
> L__u2:=diff(L,u2):
> L__u1t:=&dt L__u1:
> L__u2t:=&dt L__u2:
> eq1:=L__u1t-L__q1=0:
> eq1:=subs(q1t=u1,q2t=u2,eq1):
> eq2:=L__u2t-L__q2=0:
> eq2:=subs(q1t=u1,q2t=u2,eq2):
> kde:={q1t=u1,q2t=u2}:
> eq1:={u1t=solve(eq1,u1t)}:
> eq2:={u2t=solve(eq2,u2t)}:
> eqs:= eq1 union eq2 union kde:
> # ––––––– IV. Solver ––––––– #
> eqst:=subs(toTimeFunction,eqs):
> eqst:=subs(param,eqst):
> ff:=dsolve(eqst union InitCond,{q1(t),q2(t),u1(t),u2(t)},type=numeric,
maxfun=500000):
> #plotpath_2d(ff,r__A,param,0,10,800,black,0,N):
> #animatepath_2d(ff,r__A,param,0,10,1000,200,red,2,N):
> gg:=dsolve(eqst union InitCond,{q1(t),q2(t),u1(t),u2(t)},type=numeric,
maxfun=500000,output=Array([seq((1/4000*i),i=0..10000)])):
> gg:=gg(2,1):
> Export(gg,"2_body.xls","2 body"):
Appendix B
Matlab code
In this Appendix the Matlab code used can be found.
B.1 Planet in Keplerian orbit

B.1.1 twobody.m
clear all, close all, clc

% Read input data from Maple-generated .xls - file
A = xlsread('2_body.xls');
% Defining cartesian coordinates from polar coordinates

x = A(:,2).*cos(A(:,3));
y = A(:,2).*sin(A(:,3));
% ------------------------ %
% ---- Phase Portrait ---- %
% ------------------------ %
figure
plot(A(:,2),A(:,4),'Color',[0 0.447 0.8])
title('Phase Portrait: q_1(t)','FontSize',20)
xlabel('$$ q_1(t) $$','Interpreter','LaTex','FontSize',16)
ylabel('$$ \dot{q}_1(t) $$','Interpreter','LaTex','FontSize',16)
grid on
figure
plot(A(:,3),A(:,5),'Color',[0 0.447 0.8])
title('Phase Portrait: q_2(t)','FontSize',20)
xlabel('$$ q_2(t) $$','Interpreter','LaTex','FontSize',16)
ylabel('$$ \dot{q}_2(t) $$','Interpreter','LaTex','FontSize',16)
grid on
67
% ----------------------------- %
% ----- Kepler, First law ----- %
% ----------------------------- %
% Ellipse dimensions
a = (max(x)+min(x))/2-min(x);
b = (max(y)-min(y))/2;
e = sqrt(1-(b/a)^2);
h = (max(x)-a);
k = (max(y)-b);
areal_velocity = abs(A(:,2).^2.*A(:,5))/2; % Calculating areal velocity
% ------ Insertion of coordinates in ellipse equation ------ %

ellipse = ((x-(max(x)-a))./a).^2+((y-(max(y)-b))./b).^2;
ellipse2 = ((x-(min(x)+a))./a).^2+((y-(min(y)+b))./b).^2;
ellipse_analysis = [min(ellipse)
max(ellipse)
mean(ellipse)
median(ellipse)]; % [Min, max, mean, median]-values
ellipse_analysis2 = [min(ellipse2)
max(ellipse2)
mean(ellipse2)
median(ellipse2)]; % [Min, max, mean, median]-values
% ---------------------------- %
% ---- Kepler, Second law ---- %
% ---------------------------- %
% Physical parameters
M = 250;
m = 1;
G = 5e-1;
mu = G*(M+m);
% Positional vector, eccentric anomaly and time function (Kepler Eq.)

r = @(theta) a.*(1-e.^2)./(1+e.*cos(theta)); % Radius r(theta)
rad = @(theta) 0.5*(a.*(1-e.^2)./(1+e.*cos(theta))).^2; % (r^2)/2
E = @(theta) 2*atan(sqrt((1-e)/(1+e))*tan(theta/2)); % Eccentric anomaly
tp = @(theta) sqrt(a^3/mu)*(eccentric_anomaly_v2(theta,e)...
-e*sin(eccentric_anomaly_v2(theta,e))); % Kepler's equation
% --------------------------- %
% ---- Various variables ---- %
% --------------------------- %
% Radius for both bodies, center loc. for planet (Used in plot)
radius = 0.035;
radius2 = 0.075;
center = [x y];
theta = linspace(0,2*pi);
% Find index of one complete orbit
tmp = abs(A(:,3)-2*pi);
[idx idx] = min(tmp);
% -------------------------- %
% ---- Plot environment ---- %
% -------------------------- %
figure
% Plot the trajectory of the orbit
trajectory = plot(x,y,'k');
hold on
% Legend: Planet
X = -2.2 + radius.*cos(theta);
Y = -1.2 + radius.*sin(theta);
legend_planet=fill(X,Y,[255 0 0]/255);
text(-2.13, -1.2,'- Orbiting Planet')
% Legend: Sun
X = -0.6 + radius2.*cos(theta);
Y = -1.2 + radius2.*sin(theta);
legend_sun=fill(X,Y,[255 255 0]/255);
set(legend_sun,'EdgeColor','r');
text(-0.48, -1.2,'- Sun')
% Semi-major/minor axis
line_major = line([min(x) max(x)],[0 0],'Color',[160 160 160]/255,...
'LineWidth',1.3,'LineStyle','--');
line_minor = line([x(find(y==min(y))) x(find(y==max(y)))],[min(y) max(y)],...
'Color',[160 160 160]/255,'LineWidth',1.3,'LineStyle','--');
% Creating stationary Sun

X2 = radius2.*cos(theta);
Y2 = radius2.*sin(theta);
sun=fill(X2,Y2,[255 255 0]/255);
set(sun,'EdgeColor','r');
% ------------------------- %
% ---- Animation-block ---- %
% ------------------------- %
% Null variables/strings
area = 0; area2 = 0; j = 0; k = 0; l = 0;
time1 = '0'; time2 = '0'; time_orbit = '0';
% Iterate through all dicrete time steps.

for i = 1:length(A)
% Define axis for figure
axis([-2.8 0.5 -1.65 1.65])
% Draw circles
X = center(i,1)+radius.*cos(theta); % Create X values for the circle
Y = center(i,2)+radius.*sin(theta);
planet=fill(X,Y,[255 0 0]/255);
set(planet,'EdgeColor','None');
% Coordinates for the radius-line.
XR = linspace(0,r(A(i,3))*cos(A(i,3)));
YR = linspace(0,r(A(i,3))*sin(A(i,3)));
% Shaded area plot-environment.

if r(A(i,3))*sin(A(i,3)) >= 0
delete(time1)
if area <= pi / 12
x_draw = [0 x(1:i)'];
y_draw = [0 y(1:i)'];
A1 = patch(x_draw,y_draw,[224 224 224]/255);
uistack(A1,'bottom')
j = j + 1;
save_area_1 = area;
else
str = '$$ A_{1}=\frac{\pi}{12} $$';
text(-0.72, 0.625,str,'Interpreter','latex','fontsize',9)
end
area = integral(rad,0,A(i,3));
time1 = text(-0.8, 0.1,['{\Delta}t_{1}= ',num2str(A(j,1),...
'%.4f'),' s'],'fontsize',8);
l = l + 1;
elseif r(A(i,3))*sin(A(i,3)) < 0
delete(time2)
if area2 <= pi / 12
x_draw = [radius2*cos(A(l,3))' x(l:i)' radius2*cos(A(i,3))'];
y_draw = [radius2*sin(A(l,3))' y(l:i)' radius2*sin(A(i,3))'];
A2 = patch(x_draw,y_draw,[224 224 224]/255);
k = k + 1;
save_area_2 = area2;
uistack(A2,'bottom')
else
str = '$$ A_{2}=\frac{\pi}{12} $$';
text(-2.43, -0.32,str,'Interpreter','latex','fontsize',9)
end
area2 = integral(rad,pi,A(i,3));
time2 = text(-2, 0.1,['{\Delta}t_{2}= ',num2str(A(k+l,1)-A(l,1),...
'%.4f'),' s'],'fontsize',8);
else
end
% Only plot radius-line for theta < 2*pi
if A(i,3) < 2*pi
radius_line = line(XR,YR);
end
% Clock for elapsed time

elapsed_time = text(-1.74, 1.25,['Elapsed time: ' num2str(A(i,1),...
'%.4f') ' s'],'fontsize',12);
if A(i,3) >= 2*pi
time_orbit = text(-1.7, 1.1,['Orbital time: \tau = ' ...
num2str(A(idx,1),'%.4f') ' s'],'fontsize',10);
end
drawnow
% Delete objects to increase simulation rate.

if i < length(A)
delete(planet)
if area < pi/12
delete(A1)
end
if A(i,3) > pi && area2 < pi/12
delete(A2)
end
delete(radius_line)
delete(elapsed_time)
if A(i,3) >= 2*pi
delete(time_orbit)
end
end
end
B.1.2 areal_velocity.m
clear all, close all, clc
% Read input data from maple-generated .xls - file

A = xlsread('2_body.xls');
% Defining cartesian coordinates from polar coordinates

x = A(:,2).*cos(A(:,3));
y = A(:,2).*sin(A(:,3));
% ------------------- %
% ---- Kepler II ---- %
% ------------------- %
% Ellipse dimensions
a = (max(x)+min(x))/2-min(x);
b = max(y);
e = sqrt(1-(b/a)^2);
% Physical parameters
M = 250;
m = 1;
G = 5e-1;
mu = G*(M+m);
arealvelocity = abs(A(:,2).^2.*A(:,5))/2;
rad = @(theta) 0.5*(a.*(1-e.^2)./(1+e.*cos(theta))).^2;

tp = @(theta) sqrt(a^3/mu)*(eccentric_anomaly_v2(theta,e)...
-e*sin(eccentric_anomaly_v2(theta,e))); % Keplers equation
theta = 0:0.01:2*pi;
for i = 1:length(theta)
Area(i) = integral(rad,0,theta(i));
time(i) = tp(theta(i));
end
P = polyfit(time,Area,1);
figure
plot(time,Area,'b',time,polyval(P,time)+0.3,'r-','LineWidth',1.5)
title('Swept area as a function of elapsed time')
legend(['Calculated data'],['Reference line: A = ',num2str(P(1),'%.2f'),'t + 0.3'])
grid on
xlabel('Elapsed time since periapsis [s]')
ylabel('Swept area [area unit]')
for i = 1:length(theta)-1
Area_vel(i) = (Area(i+1)-Area(i))/(time(i+1)-time(i));
end
figure
plot(time(1:end-1),Area_vel,'b','LineWidth',1.5)
axis([0 1 4.3433 4.344])
grid on
title('Areal velocity as a function of elapsed time since periapsis')
ylabel('Areal velocity [(area unit)/s]')
figure
plot(A(1:4001,1),arealvelocity(1:4001),'b','LineWidth',1.5)
axis([0 1 4.33 4.34])
grid on
title('Areal velocity as a function of elapsed time since periapsis')
ylabel('Areal velocity [(area unit)/s]')
74 APPENDIX B. MATLAB CODE
B.1.3 eccentric_anomaly.m
% Calculates the eccentric anomaly

function E = eccentric_anomaly_v2(theta,e)
if theta > pi && theta < 2*pi

E = 2*pi + 2*atan(sqrt((1-e)/(1+e))*tan(theta/2));
else
E = 2*atan(sqrt((1-e)/(1+e))*tan(theta/2));
end

Fulltext01 PDF

Uploaded by

Copyright:

Available Formats

Fulltext01 PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Fulltext01 PDF

Uploaded by

Copyright:

Available Formats

Modelling and Simulation of Conservative

Dynamical Systems by Computer Algebra

SA115X Degree Project in Vehicle Engineering, First Level

FREDRIK DEHLIN and JACOB ASK OLSSON

Supervisor: Nicholas Apazidis

May 22, 2017

Conservative dynamical systems is modelled with Lagrangian mechanics using

See first, think later, then test.

3.4 Double Pendulum in 3D . . . . . . . . . . . . . . . . . . . . . . . . . 30

1.3 Problem Formulations

1.3.1 Double Pendulum in 2D

Figure 1.1: Rendered double pendulum in 2D.

The definition of the problem as formulated by Apazidis [1, p. 27]:

Consider a double pendulum consisting of two bars OA and AB of masses m1 and

The sensitivity analysis will be performed by comparing trajectory plots as well as

1.3.2 Double Pendulum in 3D

Figure 1.2: Rendered double pendulum in 3D.

The definition of the problem as formulated by Apazidis [1, p. 26]:

1.3.3 Planet in Keplerian Orbit

Figure 1.3: Rendered view of a planet in an elliptical orbit around a star.

2.1 Newtonian Mechanics

2.1.1 Newton’s Laws of Motion

I. ’Every body continues in its state of rest or of uniform rectilinear motion,

2.1.2 Energy, Work and Conservative Fields

2.1.3 Moment of Inertia

In the case of an object rotating around an arbitrary axis Q it can be shown

Parallel axis theorem

Ixy = Ix0 y0 − mxG yG , (2.15)

2.1.4 Coordinate Transforms

complete understanding of the particles motion it might be necessary to study its

Rotation of coordinate systems

T1 (α) = 0 cos(α) − sin(α)

T3 (γ) = − sin(γ) cos(γ) 0 .

in the new coordinate system S 0 .

2.2 Functionals and Variational Calculus

A fundamental problem in variational calculus consists of finding a real function

as a small variation in the argument of the functional, where η and ϕ is considered

2.3 Lagrangian Mechanics

The constraints of this system can be expressed as x2 + y 2 + z 2 = R2 , in the depen-

r1 , r2 , ..., rn → q1 , q2 , ..., qf , f = 1, 2, ..., 3n − Λ ’ , (2.24)

By specifying the configuration of a system by the use of generalized coordinates

The coordinate path is essentially a mapping from time to tuples2 of generalized

2.3.1 Lagrange’s Equations

L(qk , q˙k , t) = T (qk , q˙k ) − V (qk , t) , (2.29)

2.3.2 Hamilton’s Variational Principle

2.3.3 Lagrange’s Equations from Hamilton’s Principle

2.4 Stability of Dynamical Systems

2.5 Orbital Mechanics

2.5.1 Kepler’s Laws

’Kepler’s Laws are:

I. The orbit of each planet is an ellipse with the sun at a focus.

Figure 2.5: Keplerian motion of a planet around a star.

2.5.2 Mathematics of Ellipses

a2 (1 − e2 )(r2 cos2 (θ) + 2aer cos(θ) + a2 e2 ) + a2 r2 (1 − cos2 (θ)) = a4 (1 − e2 ) , (2.39)

which simplifies into

and by changing into polar coordinates

is important to remember that the x component is expressed as x = r(θ) cos(θ) − c

Figure 2.7: Defining the eccentric- and true anomaly.

Inserting (2.50) and (2.52) into (2.48)