c and ee 103 applied numerical computing and modeling note that it is

Question

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