C&EE 103

Applied Numerical Computing and Modeling

Final Project The Elastic Pendulum

www

Figure 1: Swinging spheres that are attached to springs.

Summer 2023

Introduction. This final project deals with the dynamics of an elastic pendulum, as depicted in Fig. 1. The

theoretical background of this engineering problem demonstrates various important aspects of engineer-

ing. First and foremost, it is a non-linear dynamics problem in which the system of a sphere and a spring

undergo a motion that is not restricted to small angular positions of the pendulum. Dynamics problems

play an important role in various fields of engineering, such as earthquake engineering, aerospace engineer-

ing, robotics, and many others. A basic understanding of this type of problems is therefore of the utmost

importance for many engineers. Moreover, with this project, you can also get a basic understanding of

material modeling since the spring is assumed to exhibit elastic material behavior, which is characterized

by non-dissipating energies. This makes the modeling considerably more straightforward than using, e.g.,

plastic materials. In this project, however, the material behavior is assumed to be non-linear, i.e. the force-

displacement curve does not obey Hooke's law as most materials exhibit, at least for large deformations,

essentially non-linear material behavior. It is demonstrated how this dynamics problem can be simplified

and turned into a physical-mathematical model expressed by an initial value problem. Since the true solu-

tion to this problem cannot be found, the problem will be solved numerically by the methods learned in this

course, such as multi-dimensional rootfinding methods, interpolation/approximation methods, numerical in-

tegration schemes, efficient linear systems of equations solvers, and numerical methods to solve initial value

problems. With this final project, you can therefore demonstrate your knowledge on numerical analysis to

find an approximate solution to this non-linear problem.

(Fig. 1 has been designed using an image from stock.adobe.com)

1

Theoretical framework. A general strategy to obtain a feasible engineering problem stated in terms of

a physical-mathematical model is to simplify the problem at hand following the engineering rationale "as

simple as possible, as complex as necessary"./nC&EE 103

Applied Numerical Computing and Modeling

mmmet.

L+u(t)

Figure 2: The elastic pendulum.

m

To this end, any type of friction is neglected in this project since their effects on the solution to this problem

are usually small. In more accurate models, friction, such as aerodynamical drag, needs to be taken into

account. As a consequence, the energy of the elastic pendulum considered in this project is a conserved

quantity. For the sake of simplicity, the spring is assumed to be considerably lighter than the sphere attached

to it, and therefore, the mass of the spring can be disregarded. The sphere, on the other hand, is not massless,

but it is restricted to the case that the forces acting on the sphere are assumed to be sufficiently small and

that the strength of the material is assumed to be sufficiently large so that any deformations of the sphere are

negligible. other words, is assumed that the sphere can be modeled as a rigid body. Many pendulum

problems, however, simplify the sphere even further to a point mass, and therefore, the model used in

this project allows for a more realistic modeling of the pendulum. Another important simplification deals

with the dimension of the problem. This generally three-dimensional problem is considered in two space

dimensions, and thus, any motion of the pendulum out of the plane is disregarded.

Based on these simplifying assumptions, the physical-mathematical model of the elastic pendulum can be

derived in an abstract way by considering the energies that appear in the system, which are the kinetic energy

and both the gravitational and elastic potential energies.

The kinetic energy of the sphere, as a rigid body in general plane motion, takes on the form

1

T() = =m(r) + }10(o)?.

V(t) = mgh +

Here, m is the mass of the sphere, ō is the velocity (vector) at its center (of mass), Ī = 2/5m r² is the mass

moment of inertia associated with the center of mass of the sphere of radius r, is its angular position,

and the dot indicates the first-order derivative with respect to the time t. According to Fig. 2, the velocity

(vector) can be expressed as 7 = {ù(t) (L + u(t))Ò (t)}T, where u(t) is the displacement of the spring of

length L such that L + u(t) is the deformed length of the spring at time t.

The potential energy of the system, as the sum of its gravitational and elastic parts, can be stated as

•u(t)

2

Summer 2023

F(u (t)) du./nC&EE 103

Applied Numerical Computing and Modeling

Here, g = 9.796 m/s² is the gravitational acceleration at UCLA, h is the height of the center of mass, which

can be expressed as −(L+u(t)) cos with the datum defined according to Fig. 2, and F(u(t)) describes the

generally non-linear relation between the displacement u(t) and the force F required to stretch/compress

the non-linear elastic spring by u(t). Here, we assume that F can be described by an at most quartic

polynomial, i.e. F(u(t)) ko aku(t)k. As a consequence, the elastic potential energy can be recast

into the form V₂ (t) = Σk k≤oªk/(k+ 1)u(t)k+¹. Since no force is required for a zero displacement, which

results in ao = 0, the elastic potential energy reduces to Ve(t) = Σk k=1 ak/(k + 1)u(t)k+¹ with unknown

coefficients at that need to be determined.

=

4

k=0

The equations of motion of the elastic pendulum can be derived from the so-called Lagrangian, which is

simply defined as L = T(t) – V(t). The Euler-Lagrange differential equations then follow from the two

conditional equations

d ƏL IL

dt Di du

In the first case, this results after some algebraic manipulations in the second-order, non-linear ordinary

differential equation

Y(t) =

u(t)

0(t)

= 0 and

v(t)

w (t)

ü(t) = (L + u(t))ġ(t)² + g cos 0 (t)

respectively.

d as

dt a

aku(t)k.

In the second case, this leads, again after some algebraic manipulations, to the second-order, non-linear

ordinary differential equation

‚ ƒ (Y(t)) =

2ů (t) (t) + g sin 0 (t)

(L+u(t))² + ²r²

ƏL

20

m

Ô (t)

To be able to solve these differential equations numerically with the methods learned in this course, the

second-order differential equations need to be turned into first-order differential equations. Upon introduc-

ing the variables v(t) = u(t) and w(t) = ġ(t), where v(t) is the speed of the sphere in the direction of

the spring and w(t) is the angular velocity of the sphere, the two coupled second-order, non-linear ordinary

differential equations derived above can indeed be turned into a system of four coupled first-order, non-

linear ordinary differential equations. The resulting initial value problem can then be stated in vector form

as follows: find the solution Y (t) that satisfies

4

3

k=1

Here, to and tf are the start and end times of the time interval of interest, respectively. Furthermore, the

vectors Y (t), f (Y (t)), and Yo are defined as

Summer 2023

= 0.

Y(t) = f(Y(t)) in (to, tf]

Y(to) = Yo.

-(L+u(t)).

v(t)

w (t)

(L + u(t))w(t)² + g cos 0 (t) — 1/1 Σ -₁ aku(t)k |

2v(t)w (t)+g sin(t) (L + u(t))

(L+u(t))² +r²

and Yo=

uo

00

VO

WO/nC&EE 103

Uj

Fi

Applied Numerical Computing and Modeling

0.84 0.98

0.05 0.17 0.22 0.28 0.35 0.46

5.8 7.3 8.5

0.4 2.2 4.6

0.51 0.59 0.63 0.69 0.75

9.1 10.7 11.1 12.3 12.8 13.5 14.2

1.12 1.25 1.33

1.47 1.64 1.82 2.11 2.32 2.64 2.92 3.27 3.65 4.31

14.9 15.2 15.8 16.4 16.3 16.7 17.1 17.5 18.4 18.9 20.7 21.9

24.1

Table 1: Data points (ui, Fi) measured in m and N.

Ui

Fi

Summer 2023

Problem statement. Write a MATLAB script to solve the given initial value problem numerically in the

given time interval with to = 0s and tf = 6s for the displacement u (t) and the speed v(t) of the sphere

relative to the undeformed spring and for the angular position (t) and the angular velocity w(t) of the

sphere. In other words, find the path that the sphere takes based on the data of the problem, which is given

by the undeformed length L = 1.7m of the spring, the mass m = 0.6kg of the sphere and its radius

r = 0.3 m, the gravitational acceleration g = 9.796 m/s² at UCLA, and the initial conditions up = 0m,

0o = 175°, vo= 0m/s, and wo = Orad/s. To model the non-linear elastic material behavior of the spring,

the polynomial function F(u(t)) needs to be determined based on the set of data points (ui, Fi) that can be

found in Table 1. Employ numerical methods learned in this course to find and visualize an approximate

solution to the given initial value problem for two cases. In the first case, consider one elastic pendulum,

and in the second case, consider three identical elastic pendulums that start to swing at the same time.

Part I. Consider the case in which one elastic pendulum swings based on the given initial conditions. To be

more precise, do the following:

a) To alleviate measurement errors, construct a least-squares fitting function that expresses the force

F(u) as a function of the displacement u. As mentioned earlier, take into account that a zero dis-

placement requires a zero force, and thus F(0) = ON shall be satisfied exactly. The least-squares

fitting function F(u) shall be constructed as a polynomial of degree ≤ 4. Your MATLAB script should

therefore be flexible enough to use any polynomial order between 1 and 4. Since the least-squares

method is not used for the numerical method later in this project, it is not required to derive F(u) by

means of basis functions. Instead, use a MATLAB built-in function to solve the resulting linear sys-

tem of equations Ma = b for the coefficients a that are required on the right-hand side of the given

initial value problem. Use a fourth-order polynomial to generate the figures for your report. However,

use also some of the lower degrees to investigate how the fitting function F(u) affects the numerical

results, and include a brief statement of your observations in your report.

b) Solve the given initial value problem numerically by the backward Euler method, and employ New-

ton's method with a tolerance & = 10-12 in each subinterval to solve the resulting non-linear system

of equations and to create the four sets of data points (tn, Un), (tn, On), (tn, Un), and (tn, wn) for all

n. Do so by employing two different distributions of nodes tn. The first solution shall be obtained

by a coarse distribution of nodes with as few nodes as possible. More precisely, to find the number

of subintervals ns1 of equal length h associated with the coarse distribution of nodes, use a value

of ns1 = 50 to start with, and, if necessary, increase this value successively in steps of 50 until a

numerical solution is found. The second solution shall be obtained by a fine distribution of nodes, for

which ns1 = 2000 subintervals of equal length h are used.

c) Since friction is not considered in this problem, the energy E(t) = T(t) + V(t) obtained by the

4/nC&EE 103

Applied Numerical Computing and Modeling

true but unknown solution Y (t) is a conserved quantity. In the case of the numerical solution y(t),

however, it cannot be expected that the energy E(t) is conserved, i.e. it cannot be expected that

E(t) = const. To analyze the quality of the numerical solution y(t) by means of the energy of the

numerical solution, generate the set of data points (tn, En) for all n.

Summer 2023

d) Once the numerical solution y(t) to the initial value problem is found in terms of the associated data

points, create an animation of the motion of the elastic pendulum on screen, i.e. without saving it into

a file. Compare the animations generated by the two different values of ns1 introduced in Part Ib),

and include a brief statement of your observations in your report.

e) To be able to plot the numerical results found in Parts Ib) and Ic) and given in terms of approximate

solutions of u (t), 0(t), v(t), w(t), and E(t), construct for each of the two values of n s7 the associated

piecewise quadratic Lagrangean interpolant q(t) that interpolates the data points (tn, Un), (tn, On),

(tn, Un), (tn, wn), and (tn, En) found in Parts Ib) and Ic). Construct each of the five interpolants q(t)

in terms of the quadratic Lagrangean interpolating polynomial in each of the ns1/2 subintervals for

both values of ns1.

f) Plot all data points (ui, F;) given in Table 1 and the least-squares fitting function F(u) found in Part

Ia) in one figure. Furthermore, create three more figures for each of the two values of ns1. In the

first of the three figures, plot the interpolant found in Part Ie) associated with the approximation of

E (t). In the second of the three figures, plot the interpolants found in Part Ie) associated with the

approximations of u (t) and v(t). Lastly, in the third of the three figures, plot the interpolants found

in Part Ie) associated with the approximations of (t) and w(t). This results in seven figures in total.

Compare the figures, and include a brief statement of your observations in your report.

Part II. Consider the case in which three identical elastic pendulums swing at the same time but based on

different initial conditions. To be more precise, do the following:

a) Expand your MATLAB script to include two additional elastic pendulums of the same type as the first

one in your script, i.e. all assumptions and the data of the problems are the same as in Part I except for

the initial conditions o. Choose 00 = 176° for the first additional pendulum and 0o = 177° for the

second additional pendulum. For a longer animation, expand the time interval by redefining the end

time as tf = 60 s.

b) Solve the three given initial value problems numerically by the trapezoidal method, and employ the

predictor-corrector method with a tolerance & = 10-12 in each subinterval to solve the resulting non-

linear system of equations and to create the four sets of data points (tn, Un), (tn, On), (tn, Vn), and

(tn, wn) for all n and for each of the three solutions. Do so by employing ns1 = 20000 subintervals

of equal length h for each solution, which is equivalent to ten times the fine distribution of nodes t

used in Part Ib).

c) Once all three numerical solutions y(t) to the initial value problems are found in terms of the asso-

ciated data points, create an animation of the motion of the three elastic pendulums on screen, i.e.

without saving it into a file. Compare the animations generated by the two different values of n sı,

and include a brief statement of your observations in your report.

5/nC&EE 103

Applied Numerical Computing and Modeling

Note that it is not required in this part to create any interpolant. As a consequence, it is not required to create

any plots.

Summer 2023

Instructions. This project is designed to build on the work you have done in this course, and thus, it focuses

on the numerical methods that you have learned in this course. Even though the problem is an engineering

problem, it is not required that you have a thorough understanding of all details of the modeling aspect,

particularly if you are not an engineering student. It is more important that you know how to solve the

problem numerically. Moreover, you are expected to work on the project independently, and thus, any

questions of clarification should be directed to the instructor.

Your work should consist of both a typed scientific report that fully documents the problem, your algorithmic

setup, and your solution and a MATLAB script that will fully solve the problem.

Your report should be a pdf-document that includes sections describing the basic problem setup and identi-

fication of the problem type you are solving, the methods of solution you are using (with details of how it is

implemented), results showing the numerical solution(s) via plot(s), a discussion of relevant implications of

your solution method (e.g. accuracy, stability, etc.), and your MATLAB script(s) in the appendix.

You are supposed to write your own MATLAB script (based on a template that provides a general struc-

ture of the script and some code snippets related to parts of the problem that are not related to the topics

of this course). Where possible, you may leverage the MATLAB scripts that have been provided during

the course. Some of these MATLAB scripts will be able to be used directly, others may require modifica-

tion/generalization for usage in your project. Your script should be capable of solving any problem of the

type you are given (with only slight adjustments). Do not use any built-in MATLAB functions as part of your

script (except for the visualization of your results and further smaller exceptions).

Turn in both your report and your MATLAB script for this project. Remember to use your name in the

report, the MATLAB filename, and in the MATLAB script, and upload your report and your MATLAB script

to gradescope.

6