Computing and Numerical Methods 2

Coursework

Academic Year 2020/2021

Deadline: 29th January 2021, 11pm GMT

Instructions

This coursework is divided into two parts. You are required to follow the instructions for each part detailed

below.

Part I:

This part consists of four exercises and involves programming in Matlab. You must answer each question

using the word template provided. Ensure that all figures are clear, labelled correctly and that you strictly

adhere to any word/space limits stated in the template. Do not alter the template in any way: textboxes

should not be changed and the written responses must be in “Calibri”, fontsize 12 (handwritten responses

are not accepted). Figures should be completely contained inside the textboxes. If the template calls only

for a plot, only a plot should be provided (anything else will not be marked).

Having completed the tutorial sheets, you will already have the core parts of the Matlab code required to

complete Part I. While marks will be awarded for each question, questions requring analysis carry a heavier

weight. Any discussions/analysis can be made solely based on the course material and observations made

regarding the questions within this coursework. Make sure that your answers to these parts are clear,

concise and to-the-point.

Part II:

This part requires you to write a program in C++. While the structure of the program is not prescribed,

marks will be awarded for well-structured code which is uses good programming practices and is clear to

read. You should ensure your program adheres to any requirements on the input and output format of

parameters and generated data.

The short report should be written in size 12 font. It should not exceed two A4 pages overall; any material

beyond the 2-page limit will not be marked. Make sure that your answers to these parts are clear, concise

and to-the-point. Plots should be legible and any text on figures should be no smaller than size 12 font.

Notation: Given a vector v, its derivative with respect to an independent variable t is denoted by v˙. Namely

dv

dt

= ˙v .

1

Figure 1: Diagram of the double pendulum. The diagram is not drawn to scale.

In this assignment you are going to numerically solve a system of ordinary differential equations which model the

motion of a simple double pendulum shown in Figure 1. The system consists of two weights of arbitrary mass, m1

and m2, connected by mass-less rigid rods of prescribed lengths, l1 and l2. As shown in Figure 1, the mass m1 is

considered as a point mass located at position (x1, y1) with a deflection angle θ1, whereas the second mass m2 is

located at position (x2, y2) with a deflection angle θ2. The energy of the system is given by

L = Linear Kinetic Energy + Potential Energy .

is the gravitational acceleration. Assuming conservation of energy, the system must satisfy

the Euler-Lagrange equation and can be described by the second-order differential equations (in the variables θ1

and θ2)

(m1 + m2)l1

¨θ1 + m2l2 cos(θ1 ? θ2)

¨θ2 + m2l2 sin(θ1 ? θ2)

˙θ

2

2 + (m1 + m2)g sin(θ1) = 0 ,

l2

¨θ2 + l1 cos(θ1 ? θ2)

¨θ1 ? l1 sin(θ1 ? θ2)

˙θ

2

1 + g sin(θ2) = 0 .

(1)

The (x, y)-coordinates of each pendulum can readily be expressed in terms of the lengths l1 and l2, and the angles

θ1 and θ2 as follows

x1 = l1 sin(θ1),

y1 = ?l1 cos(θ1),

x2 = l1 sin(θ1) + l2 sin(θ2),

y2 = ?l1 cos(θ1) ? l2 cos(θ2).

The initial position of the pendulum is such that θ1(0) = θ1,0 and θ2(0) = θ2,0 (the values of which will be specified

later) and the pendulum is initially at rest, i.e. ˙θ1(0) = ˙θ2(0) = 0.

2

m1 m2 l1 l2 θ1,0 θ2,0

1 1 0.01 1 0 π/4

Table 1: Parameters defining the pendulum considered in Q2).

m1 m2 l1 l2

1 1 1 1

Table 2: Parameters defining the pendulum considered in Q3)-Q4).

Part I - Numerical Methods

Q1) The two-dimensional, second order differential equation (1) can be expressed as a first order differential

equation of a larger dimension.

(a) Write the system of equations (1) in the form

M(θ1, θ2)

¨θ1

¨θ2

= g(θ1, θ2,

˙θ1,

˙θ2),

where the matrix M and vector g are smooth functions to be determined.

(b) Using the result of part (a) of Q1), express (1) as a system of first-order differential equations in terms of

a vector v of suitable dimension. Namely, express the second order differential equation (1) in the form

v˙ = f(v),

where f is a smooth function to be determined. [5%]

Q2) For this part of the problem, suppose that the various parameters characterising the pendulum are as speci-

fied in Table 1.

(a) For the initial conditions given in Table 1, solve the differential equation obtained in Q1)(b) using explicit

Euler with step size h = 0.1 for the interval t = [0, 10].

i. Plot the angles θ1 and θ2 obtained numerically as a function of time. Plot the two angles on different

axes (e.g. using two subplots).

ii. In a new figure, plot the position (x1, y1) for each time (such a plot is called a phase portrait) in red.

On the same figure, plot the position (x2, y2) for each time in blue.

iii. Inspecting the plots obtained in the preceding two sub-problems, do you believe the numerical

solution obtained represents well the behaviour of the double pendulum? Provide a brief explanation

for your answer. [15%]

Q3) For this part of the problem, suppose that the various parameters characterising the pendulum and are as

specified in Table 2.

(a) Let θ1,0 = π and θ2,0 =

π

2

. Solve the differential equation obtained in Q1)(b) using 4th order Runge-Kutta

with step size h = 0.001 for the interval t = [0, 10].

i. Plot the angles θ1 and θ2 obtained numerically as a function of time. Plot the two angles on different

axes (e.g. using two subplots).

ii. Plot x1, y1, x2 and y2 obtained numerically as a function of time. Plot each variable on different axes

(e.g. using four subplots).

iii. In a new figure, plot the position (x1, y1) for each time (such a plot is called a phase portrait) in red.

On the same figure, plot the position (x2, y2) for each time in blue.

(b) Let θ1,0 = π + 0.001 and θ2,0 =

π

2

. Solve the differential equation obtained in Q1)(b) using 4th order

Runge-Kutta with step size h = 0.001 for the interval t = [0, 10].

i. Plot the angles θ1 and θ2 obtained numerically as a function of time. Plot the two angles on different

axes (e.g. using two subplots).

3

ii. Plot x1, y1, x2 and y2 obtained numerically as a function of time. Plot each variable on different axes

(e.g. using four subplots).

iii. In a new figure, plot the position (x1, y1) for each time (such a plot is called a phase portrait) in red.

On the same figure, plot the position (x2, y2) for each time in blue.

(c) Compare the two numerical solutions obtained in Q3)(a) and Q3)(b). Comment on what you observe and

discuss briefly whether your observations are a consequence of the numerical solver used or whether

it is an artefact of the differential equation itself. [15%]

Q4) For this part of the problem, suppose that the various parameters characterising the pendulum are again as

specified in Table 2.

(a) Let θ1,0 = π and θ2,0 =

π

2

. Solve the differential equation obtained in Q1)(b) using explicit Euler with

step size h = 0.001 for the interval t = [0, 10].

axes (e.g. using two subplots).

ii. Plot x1, y1, x2 and y2 obtained numerically as a function of time. Plot each variable on different axes

(e.g. using four subplots).

iii. In a new figure, plot the position (x1, y1) for each time (such a plot is called a phase portrait) in red.

On the same figure, plot the position (x2, y2) for each time in blue.

(b) Let θ1,0 = π and θ2,0 =

π

2

. Solve the differential equation obtained in Q1) using implicit Euler with step

size h = 0.001 for the interval t = [0, 10].

axes (e.g. using two subplots).

(e.g. using four subplots).

On the same figure, plot the position (x2, y2) for each time in blue.

(c) Compare the two numerical solutions obtained in Q4)(a) and Q4)(b). Comment briefly on what you

observe and discuss whether or not one numerical solution is better than the other. [15%]

Continues on next page ...

4

Part II: C++ Programming

Q5) Write a C++ code to numerically simulate the motion of the double-pendulum. Solve the system, given by

equation (1) from time 0 to time T. Your code should:

? read a file called parameters.txt, located in the current working directory, to supply the values of m1,

m2, l1, l2, θ1,0, θ2,0, ?t, T. These values must be specified in the file in the order given above, with

values on a single line separated by spaces, for example:

1.5 2.0 1.0 1.0 1.571 2.356 0.01 10.0

? integrate the system using the 4th-order Runge-Kutta scheme.

? write the time t and (x, y)-coordinates of the masses at each time step to the file output.txt in the

current working directory;

? not require interaction with the user through the terminal during execution;

? not use third-party libraries or code beyond the standard C++ header files;

Marks will be awarded for the following:

? Code correctness;

? Code design (use of object-oriented paradigms, functions, suitable data structures);

? Code readability (layout, sensible choice of variable/function names, use of comments). [30%]

Q6) Write a short report as follows (max 2 pages).

(a) Validate your C++ code against your MATLAB implementation from Part I by generating 2D plots showing

the time evolution of the positions of the mass m2 for the following cases, given in terms of the vector

of parameters (m1, m2, l1, l2, θ1,0, θ2,0, ?t, T) as:

i. 1.0 1.0 0.01 1.0 0.0 0.785 0.0001 1.0

ii. 1.0 1.0 1.0 1.0 0.785 0.785 0.001 10.0

iii. 1.0 2.0 0.5 2.0 1.571 2.356 0.001 10.0

(b) Consider the parameter set 1.5 2.0 1.0 1.0 1.571 2.356 0.01 10.0.

i. Compare and plot the behaviour of the pendulum under this parameter set against the same parameters

but with ?t = 0.001.

ii. Identify the earliest time at which the difference between the positions of the second mass at these

different values of ?t exceeds an error threshold of = 0.05.

iii. Determine if the solution obtained using ?t = 0.001 is converged at the final time, based on the

threshold of = 0.01. Explain how you arrived at your answer.

(c) Discuss the potential advantages and disadvantages of solving this problem in C++, instead of an

interpreted language.

(d) Explain the benefits of using an object-oriented programming paradigm for this problem. [20%]

Submission

You should submit three files to Blackboard for this assignment. Make sure that these files have the format

and names as specified below.

1. Part1.pdf: Part I of this coursework must be submitted as a PDF (i.e. not as a word file).

2. Part2.pdf: The report component of Part II must be submitted a PDF (i.e. not as a word file).

3. Code.zip: The code component of Part II must be submitted as a ZIP archive. It should contain the

necessary C++ source code files only. For this, the xarchiver program, available from the Accessories

menu on the remote Linux environment can be used. Please unpack the archive file again

before submitting to check all the necessary files have been included and your code compiles.

End of assignment.

5

版权所有：留学生编程辅导网 2018 All Rights Reserved 联系方式：QQ:99515681 电子信箱：99515681@qq.com

免责声明：本站部分内容从网络整理而来，只供参考！如有版权问题可联系本站删除。