Figure 18: Simple oscillating mass.
The following text derives some of the most well-known physical problems that lead to second-order ODE models of the type addressed in this book. We consider a simple spring-mass system; thereafter extended with nonlinear spring, damping, and external excitation; a spring-mass system with sliding friction; a simple and a physical (classical) pendulum; and an elastic pendulum.
Figure 18: Simple oscillating mass.
The most fundamental mechanical vibration system is depicted in Figure 18. A body with mass \( m \) is attached to a spring and can move horizontally without friction (in the wheels). The position of the body is given by the vector \( \rpos(t) = u(t)\ii \), where \( \ii \) is a unit vector in \( x \) direction. There is only one force acting on the body: a spring force \( \F_s =-ku\ii \), where \( k \) is a constant. The point \( x=0 \), where \( u=0 \), must therefore correspond to the body's position where the spring is neither extended nor compressed, so the force vanishes.
The basic physical principle that governs the motion of the body is Newton's second law of motion: \( \F=m\acc \), where \( \F \) is the sum of forces on the body, \( m \) is its mass, and \( \acc=\ddot\rpos \) is the acceleration. We use the dot for differentiation with respect to time, which is usual in mechanics. Newton's second law simplifies here to \( -\F_s=m\ddot u\ii \), which translates to $$ -ku = m\ddot u\tp$$ Two initial conditions are needed: \( u(0)=I \), \( \dot u(0)=V \). The ODE problem is normally written as $$ \begin{equation} m\ddot u + ku = 0,\quad u(0)=I,\ \dot u(0)=V\tp \tag{119} \end{equation} $$ It is not uncommon to divide by \( m \) and introduce the frequency \( \omega = \sqrt{k/m} \): $$ \begin{equation} \ddot u + \omega^2 u = 0,\quad u(0)=I,\ \dot u(0)=V\tp \tag{120} \end{equation} $$ This is the model problem in the first part of this chapter, with the small difference that we write the time derivative of \( u \) with a dot above, while we used \( u^{\prime} \) and \( u^{\prime\prime} \) in previous parts of the book.
Since only one scalar mathematical quantity, \( u(t) \), describes the complete motion, we say that the mechanical system has one degree of freedom (DOF).
For numerical simulations it is very convenient to scale (120) and thereby get rid of the problem of finding relevant values for all the parameters \( m \), \( k \), \( I \), and \( V \). Since the amplitude of the oscillations are dictated by \( I \) and \( V \) (or more precisely, \( V/\omega \)), we scale \( u \) by \( I \) (or \( V/\omega \) if \( I=0 \)): $$ \bar u = \frac{u}{I},\quad \bar t = \frac{t}{t_c}\tp$$ The time scale \( t_c \) is normally chosen as the inverse period \( 2\pi/\omega \) or angular frequency \( 1/\omega \), most often as \( t_c=1/\omega \). Inserting the dimensionless quantities \( \bar u \) and \( \bar t \) in (120) results in the scaled problem $$ \frac{d^2\bar u}{d\bar t^2} + \bar u = 0,\quad \bar u(0)=1,\ \frac{\bar u}{\bar t}(0)=\beta = \frac{V}{I\omega},$$ where \( \beta \) is a dimensionless number. Any motion that starts from rest (\( V=0 \)) is free of parameters in the scaled model!
The typical physics of the system in Figure 18 can be described as follows. Initially, we displace the body to some position \( I \), say at rest (\( V=0 \)). After releasing the body, the spring, which is extended, will act with a force \( -kI\ii \) and pull the body to the left. This force causes an acceleration and therefore increases velocity. The body passes the point \( x=0 \), where \( u=0 \), and the spring will then be compressed and act with a force \( kx\ii \) against the motion and cause retardation. At some point, the motion stops and the velocity is zero, before the spring force \( kx\ii \) has worked long enough to push the body in positive direction. The result is that the body accelerates back and forth. As long as there is no friction forces to damp the motion, the oscillations will continue forever.
Figure 19: General oscillating system.
The mechanical system in Figure 18 can easily be extended to the more general system in Figure 19, where the body is attached to a spring and a dashpot, and also subject to an environmental force \( F(t)\ii \). The system has still only one degree of freedom since the body can only move back and forth parallel to the \( x \) axis. The spring force was linear, \( \F_s=-ku\ii \), in the section Oscillating mass attached to a spring, but in more general cases it can depend nonlinearly on the position. We therefore set \( \F_s=s(u)\ii \). The dashpot, which acts as a damper, results in a force \( \F_d \) that depends on the body's velocity \( \dot u \) and that always acts against the motion. The mathematical model of the force is written \( \F_d =f(\dot u)\ii \). A positive \( \dot u \) must result in a force acting in the positive \( x \) direction. Finally, we have the external environmental force \( \F_e = F(t)\ii \).
Newton's second law of motion now involves three forces: $$ F(t)\ii + f(\dot u)\ii - s(u)\ii = m\ddot u \ii\tp$$ The common mathematical form of the ODE problem is $$ \begin{equation} m\ddot u + f(\dot u) + s(u) = F(t),\quad u(0)=I,\ \dot u(0)=V\tp \tag{121} \end{equation} $$ This is the generalized problem treated in the last part of the present chapter, but with prime denoting the derivative instead of the dot.
The most common models for the spring and dashpot are linear: \( f(\dot u) =b\dot u \) with a constant \( b\geq 0 \), and \( s(u)=ku \) for a constant \( k \).
A specific scaling requires specific choices of \( f \), \( s \), and \( F \). Suppose we have $$ f(\dot u) = b|\dot u|\dot u,\quad s(u)=ku,\quad F(t)=A\sin(\phi t)\tp$$ We introduce dimensionless variables as usual, \( \bar u = u/u_c \) and \( \bar t = t/t_c \). The scale \( u_c \) depends both on the initial conditions and \( F \), but as time grows, the effect of the initial conditions die out and \( F \) will drive the motion. Inserting \( \bar u \) and \( \bar t \) in the ODE gives $$ m\frac{u_c}{t_c^2}\frac{d^2\bar u}{d\bar t^2} + b\frac{u_c^2}{t_c^2}\left\vert\frac{d\bar u}{d\bar t}\right\vert \frac{d\bar u}{d\bar t} + ku_c\bar u = A\sin(\phi t_c\bar t)\tp$$ We divide by \( u_c/t_c^2 \) and demand the coefficients of the \( \bar u \) and the forcing term from \( F(t) \) to have unit coefficients. This leads to the scales $$ t_c = \sqrt{\frac{m}{k}},\quad u_c = \frac{A}{k}\tp$$ The scaled ODE becomes $$ \begin{equation} \frac{d^2\bar u}{d\bar t^2} + 2\beta\left\vert\frac{d\bar u}{d\bar t}\right\vert \frac{d\bar u}{d\bar t} + \bar u = \sin(\gamma\bar t), \tag{122} \end{equation} $$ where there are two dimensionless numbers: $$ \beta = \frac{Ab}{2mk},\quad\gamma =\phi\sqrt{\frac{m}{k}}\tp$$ The \( \beta \) number measures the size of the damping term (relative to unity) and is assumed to be small, basically because \( b \) is small. The \( \phi \) number is the ratio of the time scale of free vibrations and the time scale of the forcing. The scaled initial conditions have two other dimensionless numbers as values: $$ \bar u(0) = \frac{Ik}{A},\quad \frac{d\bar u}{d\bar t}=\frac{t_c}{u_c}V = \frac{V}{A}\sqrt{mk}\tp$$
Consider a variant of the oscillating body in the section Oscillating mass attached to a spring and Figure 18: the body rests on a flat surface, and there is sliding friction between the body and the surface. Figure 20 depicts the problem.
Figure 20: Sketch of a body sliding on a surface.
The body is attached to a spring with spring force \( -s(u)\ii \).
The friction force is proportional to the normal force on the surface,
\( -mg\jj \), and given by \( -f(\dot u)\ii \), where
$$ f(\dot u) = \left\lbrace\begin{array}{ll}
-\mu mg,& \dot u < 0,\\
\mu mg, & \dot u > 0,\\
0, & \dot u=0
\end{array}\right.$$
Here, \( \mu \) is a friction coefficient. With the signum function
$$ \mbox{sign(x)} = \left\lbrace\begin{array}{ll}
-1,& x < 0,\\
1, & x > 0,\\
0, & x=0
\end{array}\right.$$
we can simply write \( f(\dot u) = \mu mg\,\hbox{sign}(\dot u) \)
(the sign function is implemented by numpy.sign).
The equation of motion becomes $$ \begin{equation} m\ddot u + \mu mg\hbox{sign}(\dot u) + s(u) = 0,\quad u(0)=I,\ \dot u(0)=V\tp \tag{123} \end{equation} $$
A washing machine is placed on four springs with efficient dampers. If the machine contains just a few clothes, the circular motion of the machine induces a sinusoidal external force and the machine will jump up and down if the frequency of the external force is close to the natural frequency of the machine and its spring-damper system.
A classical problem in mechanics is the motion of a pendulum. We first consider a simple pendulum (sometimes also called a mathematical pendulum): a small body of mass \( m \) is attached to a massless wire and can oscillate back and forth in the gravity field. Figure 21 shows a sketch of the problem.
Figure 21: Sketch of a simple pendulum.
The motion is governed by Newton's 2nd law, so we need to find expressions for the forces and the acceleration. Three forces on the body are considered: an unknown force \( S \) from the wire, the gravity force \( mg \), and an air resistance force, \( \frac{1}{2}C_D\varrho A|v|v \), hereafter called the drag force, directed against the velocity of the body. Here, \( C_D \) is a drag coefficient, \( \varrho \) is the density of air, \( A \) is the cross section area of the body, and \( v \) is the magnitude of the velocity.
We introduce a coordinate system with polar coordinates and unit vectors \( \ir \) and \( \ith \) as shown in Figure 22. The position of the center of mass of the body is $$ \rpos(t) = x_0\ii + y_0\jj + L\ir,$$ where \( \ii \) and \( \jj \) are unit vectors in the corresponding Cartesian coordinate system in the \( x \) and \( y \) directions, respectively. We have that \( \ir = \cos\theta\ii +\sin\theta\jj \).
Figure 22: Forces acting on a simple pendulum.
The forces are now expressed as follows.
The velocity of the body is found from \( \rpos \): $$ \v(t) = \dot\rpos (t) = \frac{d}{d\theta}(x_0\ii + y_0\jj + L\ir)\frac{d\theta}{dt} = L\dot\theta\ith,$$ since \( \frac{d}{d\theta}\ir = \ith \). It follows that \( v=|\v|=L\dot\theta \). The acceleration is $$ \acc(t) = \dot\v(r) = \frac{d}{dt}(L\dot\theta\ith) = L\ddot\theta\ith + L\dot\theta\frac{d\ith}{d\theta}\dot\theta = = L\ddot\theta\ith - L\dot\theta^2\ir,$$ since \( \frac{d}{d\theta}\ith = -\ir \).
Newton's 2nd law of motion becomes $$ -S\ir + mg(-\sin\theta\,\ith + \cos\theta\,\ir) - \frac{1}{2}C_D\varrho AL^2|\dot\theta|\dot\theta\,\ith = mL\ddot\theta\dot\theta\,\ith - L\dot\theta^2\ir,$$ leading to two component equations $$ \begin{align} -S + mg\cos\theta &= -L\dot\theta^2, \tag{124}\\ -mg\sin\theta - \frac{1}{2}C_D\varrho AL^2|\dot\theta|\dot\theta &= mL\ddot\theta\tp \tag{125} \end{align} $$ From (124) we get an expression for \( S=mg\cos\theta + L\dot\theta^2 \), and from (125) we get a differential equation for the angle \( \theta(t) \). This latter equation is ordered as $$ \begin{equation} m\ddot\theta + \frac{1}{2}C_D\varrho AL|\dot\theta|\dot\theta + \frac{mg}{L}\sin\theta = 0\tp \tag{126} \end{equation} $$ Two initial conditions are needed: \( \theta=\Theta \) and \( \dot\theta = \Omega \). Normally, the pendulum motion is started from rest, which means \( \Omega =0 \).
Equation (126) fits the general model used in (71) in the section Generalization: damping, nonlinearities, and excitation if we define \( u=\theta \), \( f(u^{\prime}) = \frac{1}{2}C_D\varrho AL|\dot u|\dot u \), \( s(u) = L^{-1}mg\sin u \), and \( F=0 \). If the body is a sphere with radius \( R \), we can take \( C_D=0.4 \) and \( A=\pi R^2 \). Exercise 1.25: Simulate a simple pendulum asks you to scale the equations and carry out specific simulations with this model.
The motion of a compound or physical pendulum where the wire is a rod with mass, can be modeled very similarly. The governing equation is \( I\acc = \boldsymbol{T} \) where \( I \) is the moment of inertia of the entire body about the point \( (x_0,y_0) \), and \( \boldsymbol{T} \) is the sum of moments of the forces with respect to \( (x_0,y_0) \). The vector equation reads $$ \rpos\times(-S\ir + mg(-\sin\theta\ith + \cos\theta\ir) - \frac{1}{2}C_D\varrho AL^2|\dot\theta|\dot\theta\ith) = I(L\ddot\theta\dot\theta\ith - L\dot\theta^2\ir)\tp$$ The component equation in \( \ith \) direction gives the equation of motion for \( \theta(t) \): $$ \begin{equation} I\ddot\theta + \frac{1}{2}C_D\varrho AL^3|\dot\theta|\dot\theta + mgL\sin\theta = 0\tp \tag{127} \end{equation} $$
Usually one plots the mathematical quantities as functions of time to visualize the solution of ODE models. Exercise 1.25: Simulate a simple pendulum asks you to do this for the motion of a pendulum in the previous section. However, sometimes it is more instructive to look at other types of visualizations. For example, we have the pendulum and the free body diagram in Figures 21 and 22. We may think of these figures as animations in time instead. Especially the free body diagram will show both the motion of the pendulum and the size of the forces during the motion. The present section exemplifies how to make such a dynamic body diagram.
The drag force is magnified 5 times!
Dynamic physical sketches, coupled to the numerical solution of differential equations, requires a program to produce a sketch for the situation at each time level. Pysketcher is such a tool. In fact (and not surprising!) Figures 21 and 22 were drawn using Pysketcher. The details of the drawings are explained in the Pysketcher tutorial. Here, we outline how this type of sketch can be used to create an animated free body diagram during the motion of a pendulum.
Pysketcher is actually a layer of useful abstractions on top of standard plotting packages. This means that we in fact apply Matplotlib to make the animated free body diagram, but instead of dealing with a wealth of detailed Matplotlib commands, we can express the drawing in terms of more high-level objects, e.g., objects for the wire, angle \( \theta \), body with mass \( m \), arrows for forces, etc. When the position of these objects are given through variables, we can just couple those variables to the dynamic solution of our ODE and thereby make a unique drawing for each \( \theta \) value in a simulation.
Let us start with the most familiar part of the current problem: writing the solver function. We use Odespy for this purpose. We also work with dimensionless equations. Since \( \theta \) can be viewed as dimensionless, we only need to introduce a dimensionless time, here taken as \( \bar t = t/\sqrt{L/g} \). The resulting dimensionless mathematical model for \( \theta \), the dimensionless angular velocity \( \omega \), the dimensionless wire force \( \bar S \), and the dimensionless drag force \( \bar D \) is then $$ \begin{align} \frac{d\omega}{d\bar t} &= - \alpha|\omega|\omega - \sin\theta, \tag{128}\\ \frac{d\theta}{d\bar t} &= \omega, \tag{129}\\ \bar S &= \omega^2 + \cos\theta, \tag{130}\\ \bar D &= -\alpha |\omega|\omega, \tag{131} \end{align} $$ with $$ \alpha = \frac{C_D\varrho\pi R^2L}{2m}\tp$$ as a dimensionless parameter expressing the ratio of the drag force and the gravity force. The dimensionless \( \omega \) is made non-dimensional by the time, so \( \omega\sqrt{L/g} \) is the corresponding angular frequency with dimensions.
A suitable function for computing (128)-(131) is listed below.
def simulate(alpha, Theta, dt, T):
import odespy
def f(u, t, alpha):
omega, theta = u
return [-alpha*omega*abs(omega) - sin(theta),
omega]
import numpy as np
Nt = int(round(T/float(dt)))
t = np.linspace(0, Nt*dt, Nt+1)
solver = odespy.RK4(f, f_args=[alpha])
solver.set_initial_condition([0, Theta])
u, t = solver.solve(
t, terminate=lambda u, t, n: abs(u[n,1]) < 1E-3)
omega = u[:,0]
theta = u[:,1]
S = omega**2 + np.cos(theta)
drag = -alpha*np.abs(omega)*omega
return t, theta, omega, S, drag
The sketch function below applies Pysketcher objects to build
a diagram like that in Figure 22,
except that we have removed the rotation point \( (x_0,y_0) \) and
the unit vectors in polar coordinates as these objects are not
important for an animated free body diagram.
import sys
try:
from pysketcher import *
except ImportError:
print 'Pysketcher must be installed from'
print 'https://github.com/hplgit/pysketcher'
sys.exit(1)
# Overall dimensions of sketch
H = 15.
W = 17.
drawing_tool.set_coordinate_system(
xmin=0, xmax=W, ymin=0, ymax=H,
axis=False)
def sketch(theta, S, mg, drag, t, time_level):
"""
Draw pendulum sketch with body forces at a time level
corresponding to time t. The drag force is in
drag[time_level], the force in the wire is S[time_level],
the angle is theta[time_level].
"""
import math
a = math.degrees(theta[time_level]) # angle in degrees
L = 0.4*H # Length of pendulum
P = (W/2, 0.8*H) # Fixed rotation point
mass_pt = path.geometric_features()['end']
rod = Line(P, mass_pt)
mass = Circle(center=mass_pt, radius=L/20.)
mass.set_filled_curves(color='blue')
rod_vec = rod.geometric_features()['end'] - \
rod.geometric_features()['start']
unit_rod_vec = unit_vec(rod_vec)
mass_symbol = Text('$m$', mass_pt + L/10*unit_rod_vec)
rod_start = rod.geometric_features()['start'] # Point P
vertical = Line(rod_start, rod_start + point(0,-L/3))
def set_dashed_thin_blackline(*objects):
"""Set linestyle of objects to dashed, black, width=1."""
for obj in objects:
obj.set_linestyle('dashed')
obj.set_linecolor('black')
obj.set_linewidth(1)
set_dashed_thin_blackline(vertical)
set_dashed_thin_blackline(rod)
angle = Arc_wText(r'$\theta$', rod_start, L/6, -90, a,
text_spacing=1/30.)
magnitude = 1.2*L/2 # length of a unit force in figure
force = mg[time_level] # constant (scaled eq: about 1)
force *= magnitude
mg_force = Force(mass_pt, mass_pt + force*point(0,-1),
'', text_pos='end')
force = S[time_level]
force *= magnitude
rod_force = Force(mass_pt, mass_pt - force*unit_vec(rod_vec),
'', text_pos='end',
text_spacing=(0.03, 0.01))
force = drag[time_level]
force *= magnitude
air_force = Force(mass_pt, mass_pt -
force*unit_vec((rod_vec[1], -rod_vec[0])),
'', text_pos='end',
text_spacing=(0.04,0.005))
body_diagram = Composition(
{'mg': mg_force, 'S': rod_force, 'air': air_force,
'rod': rod, 'body': mass
'vertical': vertical, 'theta': angle,})
body_diagram.draw(verbose=0)
drawing_tool.savefig('tmp_%04d.png' % time_level, crop=False)
# (No cropping: otherwise movies will be very strange!)
It now remains to couple the simulate and sketch functions.
We first run simulate:
from math import pi, radians, degrees
import numpy as np
alpha = 0.4
period = 2*pi # Use small theta approximation
T = 12*period # Simulate for 12 periods
dt = period/40 # 40 time steps per period
a = 70 # Initial amplitude in degrees
Theta = radians(a)
t, theta, omega, S, drag = simulate(alpha, Theta, dt, T)
for time_level, t_ in enumerate(t):
sketch(theta, S, mg, drag, t_, time_level)
sketch function) saved in files
with names tmp_%04d.png. These can be combined to videos using
(e.g.) ffmpeg. A complete function animate for running the
simulation and creating video files is
listed below.
def animate():
# Clean up old plot files
import os, glob
for filename in glob.glob('tmp_*.png') + glob.glob('movie.*'):
os.remove(filename)
# Solve problem
from math import pi, radians, degrees
import numpy as np
alpha = 0.4
period = 2*pi # Use small theta approximation
T = 12*period # Simulate for 12 periods
dt = period/40 # 40 time steps per period
a = 70 # Initial amplitude in degrees
Theta = radians(a)
t, theta, omega, S, drag = simulate(alpha, Theta, dt, T)
# Visualize drag force 5 times as large
drag *= 5
mg = np.ones(S.size) # Gravity force (needed in sketch)
# Draw animation
import time
for time_level, t_ in enumerate(t):
sketch(theta, S, mg, drag, t_, time_level)
time.sleep(0.2) # Pause between each frame on the screen
# Make videos
prog = 'ffmpeg'
filename = 'tmp_%04d.png'
fps = 6
codecs = {'flv': 'flv', 'mp4': 'libx264',
'webm': 'libvpx', 'ogg': 'libtheora'}
for ext in codecs:
lib = codecs[ext]
cmd = '%(prog)s -i %(filename)s -r %(fps)s ' % vars()
cmd += '-vcodec %(lib)s movie.%(ext)s' % vars()
print(cmd)
os.system(cmd)
Consider a pendulum as in Figure 21, but this time the wire is elastic. The length of the wire when it is not stretched is \( L_0 \), while \( L(t) \) is the stretched length at time \( t \) during the motion.
Stretching the elastic wire a distance \( \Delta L \) gives rise to a spring force \( k\Delta L \) in the opposite direction of the stretching. Let \( \normalvec \) be a unit normal vector along the wire from the point \( \rpos_0=(x_0,y_0) \) and in the direction of \( \ith \), see Figure 22 for definition of \( (x_0,y_0) \) and \( \ith \). Obviously, we have \( \normalvec=\ith \), but in this modeling of an elastic pendulum we do not need polar coordinates. Instead, it is more straightforward to develop the equation in Cartesian coordinates.
A mathematical expression for \( \normalvec \) is $$ \normalvec = \frac{\rpos-\rpos_0}{L(t)},$$ where \( L(t)=||\rpos-\rpos_0|| \) is the current length of the elastic wire. The position vector \( \rpos \) in Cartesian coordinates reads \( \rpos(t) = x(t)\ii + y(t)\jj \), where \( \ii \) and \( \jj \) are unit vectors in the \( x \) and \( y \) directions, respectively. It is convenient to introduce the Cartesian components \( n_x \) and \( n_y \) of the normal vector: $$ \normalvec = \frac{\rpos-\rpos_0}{L(t)} = \frac{x(t)-x_0}{L(t)}\ii + \frac{y(t)-y_0}{L(t)}\jj = n_x\ii + n_y\jj\tp$$
The stretch \( \Delta L \) in the wire is $$ \Delta t = L(t) - L_0\tp$$ The force in the wire is then \( -S\normalvec=-k\Delta L\normalvec \).
The other forces are the gravity and the air resistance, just as in Figure 22. For motion in air we can neglect the added mass and buoyancy effects. The main difference is that we have a model for \( S \) in terms of the motion (as soon as we have expressed \( \Delta L \) by \( \rpos \)). For simplicity, we drop the air resistance term (but Exercise 1.27: Simulate an elastic pendulum with air resistance asks you to include it).
Newton's second law of motion applied to the body now results in $$ \begin{equation} m\ddot\rpos = -k(L-L_0)\normalvec - mg\jj \tag{132} \end{equation} $$ The two components of (132) are $$ \begin{align} \ddot x &= -\frac{k}{m}(L-L_0)n_x, \tag{133}\\ \tag{134} \\ \ddot y &= - \frac{k}{m}(L-L_0)n_y - g \tag{135}\tp \end{align} $$
Note that the derivation of the ODEs for an elastic pendulum is more straightforward than for a classical, non-elastic pendulum, since we avoid the details with polar coordinates, but instead work with Newton's second law directly in Cartesian coordinates. The reason why we can do this is that the elastic pendulum undergoes a general two-dimensional motion where all the forces are known or expressed as functions of \( x(t) \) and \( y(t) \), such that we get two ordinary differential equations. The motion of the non-elastic pendulum, on the other hand, is constrained: the body has to move along a circular path, and the force \( S \) in the wire is unknown.
The non-elastic pendulum therefore leads to a differential-algebraic equation, i.e., ODEs for \( x(t) \) and \( y(t) \) combined with an extra constraint \( (x-x_0)^2 + (y-y_0)^2 = L^2 \) ensuring that the motion takes place along a circular path. The extra constraint (equation) is compensated by an extra unknown force \( -S\normalvec \). Differential-algebraic equations are normally hard to solve, especially with pen and paper. Fortunately, for the non-elastic pendulum we can do a trick: in polar coordinates the unknown force \( S \) appears only in the radial component of Newton's second law, while the unknown degree of freedom for describing the motion, the angle \( \theta(t) \), is completely governed by the asimuthal component. This allows us to decouple the unknowns \( S \) and \( \theta \). But this is a kind of trick and not a widely applicable method. With an elastic pendulum we use straightforward reasoning with Newton's 2nd law and arrive at a standard ODE problem that (after scaling) is easy solve on a computer.
What is the initial position of the body? We imagine that first the pendulum hangs in equilibrium in its vertical position, and then it is displaced an angle \( \Theta \). The equilibrium position is governed by the ODEs with the accelerations set to zero. The \( x \) component leads to \( x(t)=x_0 \), while the \( y \) component gives $$ 0 = - \frac{k}{m}(L-L_0)n_y - g = \frac{k}{m}(L(0)-L_0) - g\quad\Rightarrow\quad L(0) = L_0 + mg/k,$$ since \( n_y=-11 \) in this position. The corresponding \( y \) value is then from \( n_y=-1 \): $$ y(t) = y_0 - L(0) = y_0 - (L_0 + mg/k)\tp$$ Let us now choose \( (x_0,y_0) \) such that the body is at the origin in the equilibrium position: $$ x_0 =0,\quad y_0 = L_0 + mg/k\tp$$
Displacing the body an angle \( \Theta \) to the right leads to the initial position $$ x(0)=(L_0+mg/k)\sin\Theta,\quad y(0)=(L_0+mg/k)(1-\cos\Theta)\tp$$ The initial velocities can be set to zero: \( x'(0)=y'(0)=0 \).
We can summarize all the equations as follows: $$ \begin{align*} \ddot x &= -\frac{k}{m}(L-L_0)n_x, \\ \ddot y &= -\frac{k}{m}(L-L_0)n_y - g, \\ L &= \sqrt{(x-x_0)^2 + (y-y_0)^2}, \\ n_x &= \frac{x-x_0}{L}, \\ n_y &= \frac{y-y_0}{L}, \\ x(0) &= (L_0+mg/k)\sin\Theta, \\ x'(0) &= 0, \\ y(0) & =(L_0+mg/k)(1-\cos\Theta), \\ y'(0) &= 0\tp \end{align*} $$ We insert \( n_x \) and \( n_y \) in the ODEs: $$ \begin{align} \ddot x &= -\frac{k}{m}\left(1 -\frac{L_0}{L}\right)(x-x_0), \tag{136}\\ \ddot y &= -\frac{k}{m}\left(1 -\frac{L_0}{L}\right)(y-y_0) - g, \tag{137}\\ L &= \sqrt{(x-x_0)^2 + (y-y_0)^2}, \tag{138}\\ x(0) &= (L_0+mg/k)\sin\Theta, \tag{139}\\ x'(0) &= 0, \tag{140}\\ y(0) & =(L_0+mg/k)(1-\cos\Theta), \tag{141}\\ y'(0) &= 0\tp \tag{142} \end{align} $$
The elastic pendulum model can be used to study both an elastic pendulum and a classic, non-elastic pendulum. The latter problem is obtained by letting \( k\rightarrow\infty \). Unfortunately, a serious problem with the ODEs (136)-(137) is that for large \( k \), we have a very large factor \( k/m \) multiplied by a very small number \( 1-L_0/L \), since for large \( k \), \( L\approx L_0 \) (very small deformations of the wire). The product is subject to significant round-off errors for many relevant physical values of the parameters. To circumvent the problem, we introduce a scaling. This will also remove physical parameters from the problem such that we end up with only one dimensionless parameter, closely related to the elasticity of the wire. Simulations can then be done by setting just this dimensionless parameter.
The characteristic length can be taken such that in equilibrium, the scaled length is unity, i.e., the characteristic length is \( L_0+mg/k \): $$ \bar x = \frac{x}{L_0+mg/k},\quad \bar y = \frac{y}{L_0+mg/k}\tp$$ We must then also work with the scaled length \( \bar L = L/(L_0+mg/k) \).
Introducing \( \bar t=t/t_c \), where \( t_c \) is a characteristic time we have to decide upon later, one gets $$ \begin{align*} \frac{d^2\bar x}{d\bar t^2} &= -t_c^2\frac{k}{m}\left(1 -\frac{L_0}{L_0+mg/k}\frac{1}{\bar L}\right)\bar x,\\ \frac{d^2\bar y}{d\bar t^2} &= -t_c^2\frac{k}{m}\left(1 -\frac{L_0}{L_0+mg/k}\frac{1}{\bar L}\right)(\bar y-1) -t_c^2\frac{g}{L_0 + mg/k},\\ \bar L &= \sqrt{\bar x^2 + (\bar y-1)^2},\\ \bar x(0) &= \sin\Theta,\\ \bar x'(0) &= 0,\\ \bar y(0) & = 1 - \cos\Theta,\\ \bar y'(0) &= 0\tp \end{align*} $$ For a non-elastic pendulum with small angles, we know that the frequency of the oscillations are \( \omega = \sqrt{L/g} \). It is therefore natural to choose a similar expression here, either the length in the equilibrium position, $$ t_c^2 = \frac{L_0+mg/k}{g}\tp$$ or simply the unstretched length, $$ t_c^2 = \frac{L_0}{g}\tp$$ These quantities are not very different (since the elastic model is valid only for quite small elongations), so we take the latter as it is the simplest one.
The ODEs become $$ \begin{align*} \frac{d^2\bar x}{d\bar t^2} &= -\frac{L_0k}{mg}\left(1 -\frac{L_0}{L_0+mg/k}\frac{1}{\bar L}\right)\bar x,\\ \frac{d^2\bar y}{d\bar t^2} &= -\frac{L_0k}{mg}\left(1 -\frac{L_0}{L_0+mg/k}\frac{1}{\bar L}\right)(\bar y-1) -\frac{L_0}{L_0 + mg/k},\\ \bar L &= \sqrt{\bar x^2 + (\bar y-1)^2}\tp \end{align*} $$ We can now identify a dimensionless number $$ \beta = \frac{L_0}{L_0 + mg/k} = \frac{1}{1+\frac{mg}{L_0k}},$$ which is the ratio of the unstretched length and the stretched length in equilibrium. The non-elastic pendulum will have \( \beta =1 \) (\( k\rightarrow\infty \)). With \( \beta \) the ODEs read $$ \begin{align} \frac{d^2\bar x}{d\bar t^2} &= -\frac{\beta}{1-\beta}\left(1- \frac{\beta}{\bar L}\right)\bar x, \tag{143}\\ \frac{d^2\bar y}{d\bar t^2} &= -\frac{\beta}{1-\beta}\left(1- \frac{\beta}{\bar L}\right)(\bar y-1) -\beta, \tag{144}\\ \bar L &= \sqrt{\bar x^2 + (\bar y-1)^2}, \tag{145}\\ \bar x(0) &= (1+\epsilon)\sin\Theta, \tag{146}\\ \frac{d\bar x}{d\bar t}(0) &= 0, \tag{147}\\ \bar y(0) &= 1 - (1+\epsilon)\cos\Theta, \tag{148}\\ \frac{d\bar y}{d\bar t}(0) &= 0, \tag{149} \end{align} $$ We have here added a parameter \( \epsilon \), which is an additional downward stretch of the wire at \( t=0 \). This parameter makes it possible to do a desired test: vertical oscillations of the pendulum. Without \( \epsilon \), starting the motion from \( (0,0) \) with zero velocity will result in \( x=y=0 \) for all times (also a good test!), but with an initial stretch so the body's position is \( (0,\epsilon) \), we will have oscillatory vertical motion with amplitude \( \epsilon \) (see Exercise 1.26: Simulate an elastic pendulum).
We immediately see that as \( k\rightarrow\infty \) (i.e., we obtain a non-elastic pendulum), \( \beta\rightarrow 1 \), \( \bar L\rightarrow 1 \), and we have very small values \( 1-\beta\bar L^{-1} \) divided by very small values \( 1-\beta \) in the ODEs. However, it turns out that we can set \( \beta \) very close to one and obtain a path of the body that within the visual accuracy of a plot does not show any elastic oscillations. (Should the division of very small values become a problem, one can study the limit by L'Hospital's rule: $$ \lim_{\beta\rightarrow 1}\frac{1 - \beta \bar L^{-1}}{1-\beta} = \frac{1}{\bar L},$$ and use the limit \( \bar L^{-1} \) in the ODEs for \( \beta \) values very close to 1.)
Figure 23: Sketch of one-wheel vehicle on a bumpy road.
We consider a very simplistic vehicle, on one wheel, rolling along a bumpy road. The oscillatory nature of the road will induce an external forcing on the spring system in the vehicle and cause vibrations. Figure 23 outlines the situation.
To derive the equation that governs the motion, we must first establish the position vector of the black mass at the top of the spring. Suppose the spring has length \( L \) without any elongation or compression, suppose the radius of the wheel is \( R \), and suppose the height of the black mass at the top is \( H \). With the aid of the \( \rpos_0 \) vector in Figure 23, the position \( \rpos \) of the center point of the mass is $$ \begin{equation} \rpos = \rpos_0 + 2R\jj + L\jj + u\jj + \half H\jj,\ \tag{150} \end{equation} $$ where \( u \) is the elongation or compression in the spring according to the (unknown and to be computed) vertical displacement \( u \) relative to the road. If the vehicle travels with constant horizontal velocity \( v \) and \( h(x) \) is the shape of the road, then the vector \( \rpos_0 \) is $$ \rpos_0 = vt\ii + h(vt)\jj,$$ if the motion starts from \( x=0 \) at time \( t=0 \).
The forces on the mass is the gravity, the spring force, and an optional damping force that is proportional to the vertical velocity \( \dot u \). Newton's second law of motion then tells that $$ m\ddot\rpos = -mg\jj - s(u) - b\dot u\jj\tp$$ This leads to $$ m\ddot u = - s(u) - b\dot u - mg -mh''(vt)v^2 $$
To simplify a little bit, we omit the gravity force \( mg \) in comparison with the other terms. Introducing \( u' \) for \( \dot u \) then gives a standard damped, vibration equation with external forcing: $$ \begin{equation} mu'' + bu' + s(u) = -mh''(vt)v^2\tp \tag{151} \end{equation} $$ Since the road is normally known just as a set of array values, \( h'' \) must be computed by finite differences. Let \( \Delta x \) be the spacing between measured values \( h_i= h(i\Delta x) \) on the road. The discrete second-order derivative \( h'' \) reads $$ q_i = \frac{h_{i-1} - 2h_i + h_{i+1}}{\Delta x^2}, \quad i=1,\ldots,N_x-1\tp$$ We may for maximum simplicity set the end points as \( q_0=q_1 \) and \( q_{N_x}=q_{N_x-1} \). The term \( -mh''(vt)v^2 \) corresponds to a force with discrete time values $$ F^n = -mq_n v^2,\quad \Delta t = v^{-1}\Delta x\tp$$ This force can be directly used in a numerical model $$ [mD_tD_t u + bD_{2t} u + s(u) = F]^n\tp$$ Software for computing \( u \) and also making an animated sketch of the motion like we did in the section Dynamic free body diagram during pendulum motion is found in a separate project on the web: https://github.com/hplgit/bumpy. You may start looking at the tutorial.
A bouncing ball is a ball in free vertically fall until it impacts the ground, but during the impact, some kinetic energy is lost, and a new motion upwards with reduced velocity starts. After the motion is retarded, a new free fall starts, and the process is repeated. At some point the velocity close to the ground is so small that the ball is considered to be finally at rest.
The motion of the ball falling in air is governed by Newton's second law \( F=ma \), where \( a \) is the acceleration of the body, \( m \) is the mass, and \( F \) is the sum of all forces. Here, we neglect the air resistance so that gravity \( -mg \) is the only force. The height of the ball is denoted by \( h \) and \( v \) is the velocity. The relations between \( h \), \( v \), and \( a \), $$ h'(t)= v(t),\quad v'(t) = a(t),$$ combined with Newton's second law gives the ODE model $$ \begin{equation} h^{\prime\prime}(t) = -g, \tag{152} \end{equation} $$ or expressed alternatively as a system of first-order equations: $$ \begin{align} v'(t) &= -g, \tag{153} \\ h'(t) &= v(t)\tp \tag{154} \end{align} $$ These equations govern the motion as long as the ball is away from the ground by a small distance \( \epsilon_h > 0 \). When \( h < \epsilon_h \), we have two cases.
Consider two astronomical objects \( A \) and \( B \) that attract each other by gravitational forces. \( A \) and \( B \) could be two stars in a binary system, a planet orbiting a star, or a moon orbiting a planet. Each object is acted upon by the gravitational force due to the other object. Consider motion in a plane (for simplicity) and let \( (x_A,y_A) \) and \( (x_B,y_B) \) be the positions of object \( A \) and \( B \), respectively.
Newton's second law of motion applied to each object is all we need to set up a mathematical model for this physical problem: $$ \begin{align} m_A\ddot\x_A &= \F, \tag{155}\\ m_B\ddot\x_B &= -\F, \tag{156} \end{align} $$ where \( F \) is the gravitational force $$ \F = \frac{Gm_Am_B}{||\rpos||^3}\rpos,$$ where $$ \rpos(t) = \x_B(t) - \x_A(t),$$ and \( G \) is the gravitational constant: \( G=6.674\cdot 10^{-11}\hbox{ Nm}^2/\hbox{kg}^2 \).
A problem with these equations is that the parameters are very large (\( m_A \), \( m_B \), \( ||\rpos|| \)) or very small (\( G \)). The rotation time for binary stars can be very small and large as well. It is therefore advantageous to scale the equations. A natural length scale could be the initial distance between the objects: \( L=\rpos(0) \). We write the dimensionless quantities as $$ \bar\x_A = \frac{\x_A}{L},\quad\bar\x_B = \frac{\x_B}{L},\quad \bar t = \frac{t}{t_c}\tp$$ The gravity force is transformed to $$ \F = \frac{Gm_Am_B}{L^2||\bar\rpos||^3}\bar\rpos,\quad \bar\rpos = \bar\x_B - \bar\x_A,$$ so the first ODE for \( \x_A \) becomes $$ \frac{d^2 \bar\x_A}{d\bar t^2} = \frac{Gm_Bt_c^2}{L^3}\frac{\bar\rpos}{||\bar\rpos||^3}\tp$$ Assuming that quantities with a bar and their derivatives are around unity in size, it is natural to choose \( t_c \) such that the fraction \( Gm_Bt_c/L^2=1 \): $$ t_c = \sqrt{\frac{L^3}{Gm_B}}\tp$$ From the other equation for \( \x_B \) we get another candidate for \( t_c \) with \( m_A \) instead of \( m_B \). Which mass we choose play a role if \( m_A\ll m_B \) or \( m_B\ll m_A \). One solution is to use the sum of the masses: $$ t_c = \sqrt{\frac{L^3}{G(m_A+m_B)}}\tp$$ Taking a look at Kepler's laws of planetary motion, the orbital period for a planet around the star is given by the \( t_c \) above, except for a missing factor of \( 2\pi \), but that means that \( t_c^{-1} \) is just the angular frequency of the motion. Our characteristic time \( t_c \) is therefore highly relevant. Introducing the dimensionless number $$ \alpha = \frac{m_A}{m_B},$$ we can write the dimensionless ODE as $$ \begin{align} \frac{d^2 \bar\x_A}{d\bar t^2} &= \frac{1}{1+\alpha}\frac{\bar\rpos}{||\bar\rpos||^3}, \tag{157}\\ \frac{d^2 \bar\x_B}{d\bar t^2} &= \frac{1}{1+\alpha^{-1}}\frac{\bar\rpos}{||\bar\rpos||^3}\tp \tag{158} \end{align} $$
In the limit \( m_A\ll m_B \), i.e., \( \alpha\ll 1 \), object B stands still, say \( \bar\x_B=0 \), and object A orbits according to $$ \frac{d^2 \bar\x_A}{d\bar t^2} = -\frac{\bar\x_A}{||\bar \x_A||^3}\tp$$
To better see the motion, and that our scaling is reasonable, we introduce polar coordinates \( r \) and \( \theta \): $$ \bar\x_A = r\cos\theta\ii + r\sin\theta\jj,$$ which means \( \bar\x_A \) can be written as \( \bar\x_A =r\ir \). Since $$ \frac{d}{dt}\ir = \dot\theta\ith,\quad \frac{d}{dt}\ith = -\dot\theta\ir,$$ we have $$ \frac{d^2 \bar\x_A}{d\bar t^2} = (\ddot r - r\dot\theta^2)\ir + (r\ddot\theta + 2\dot r\dot\theta)\ith\tp $$ The equation of motion for mass A is then $$ \begin{align*} \ddot r - r\dot\theta^2 &= -\frac{1}{r^2},\\ r\ddot\theta + 2\dot r\dot\theta &= 0\tp \end{align*} $$ The special case of circular motion, \( r=1 \), fulfills the equations, since the latter equation then gives \( \dot\theta =\hbox{const} \) and the former then gives \( \dot\theta = 1 \), i.e., the motion is \( r(t)=1 \), \( \theta(t)=t \), with unit angular frequency as expected and period \( 2\pi \) as expected.
Although the term "mechanical vibrations" is used in the present book, we must mention that the same type of equations arise when modeling electric circuits. The current \( I(t) \) in a circuit with an inductor with inductance \( L \), a capacitor with capacitance \( C \), and overall resistance \( R \), is governed by $$ \begin{equation} \ddot I + \frac{R}{L}\dot I + \frac{1}{LC}I = \dot V(t), \tag{159} \end{equation} $$ where \( V(t) \) is the voltage source powering the circuit. This equation has the same form as the general model considered in the section Generalization: damping, nonlinearities, and excitation if we set \( u=I \), \( f(u^{\prime})=bu^{\prime} \) and define \( b=R/L \), \( s(u) = L^{-1}C^{-1}u \), and \( F(t)=\dot V(t) \).
We consider the scaled ODE model (122) from the section General mechanical vibrating system. After scaling, the amplitude of \( u \) will have a size about unity as time grows and the effect of the initial conditions die out due to damping. However, as \( \gamma\rightarrow 1 \), the amplitude of \( u \) increases, especially if \( \beta \) is small. This effect is called resonance. The purpose of this exercise is to explore resonance.
a)
Figure out how the solver function in vib.py can be called
for the scaled ODE (122).
b) Run \( \gamma =5, 1.5, 1.1, 1 \) for \( \beta=0.005, 0.05, 0.2 \). For each \( \beta \) value, present an image with plots of \( u(t) \) for the four \( \gamma \) values.
Filename: resonance.
Consider a sliding box on a flat surface as modeled in the section A sliding mass attached to a spring. As spring force we choose the nonlinear formula $$ s(u) = \frac{k}{\alpha}\tanh(\alpha u) = ku + \frac{1}{3}\alpha^2 ku^3 + \frac{2}{15}\alpha^4 k u^5 + \Oof{u^6}\tp$$
a) Plot \( g(u)=\alpha^{-1}\tanh(\alpha u) \) for various values of \( \alpha \). Assume \( u\in [-1,1] \).
b) Scale the equations using \( I \) as scale for \( u \) and \( \sqrt{m/k} \) as time scale.
c) Implement the scaled model in b). Run it for some values of the dimensionless parameters.
Filename: sliding_box.
The section Bouncing ball presents a model for a bouncing ball. Choose one of the two ODE formulation, (152) or (153)-(154), and simulate the motion of a bouncing ball. Plot \( h(t) \). Think about how to plot \( v(t) \).
Hint. A naive implementation may get stuck in repeated impacts for large time step sizes. To avoid this situation, one can introduce a state variable that holds the mode of the motion: free fall, impact, or rest. Two consecutive impacts imply that the motion has stopped.
Filename: bouncing_ball.
Simulation of simple pendulum can be carried out by using the mathematical model derived in the section Motion of a pendulum and calling up functionality in the vib.py file (i.e., solve the second-order ODE by centered finite differences).
a) Scale the model. Set up the dimensionless governing equation for \( \theta \) and expressions for dimensionless drag and wire forces.
b)
Write a function for computing
\( \theta \) and the dimensionless drag force and the force in the wire,
using the solver function in
the vib.py file. Plot these three quantities
below each other (in subplots) so the graphs can be compared.
Run two cases, first one in the limit of \( \Theta \) small and
no drag, and then a second one with \( \Theta=40 \) degrees and \( \alpha=0.8 \).
Filename: simple_pendulum.
The section Motion of an elastic pendulum describes a model for an elastic pendulum, resulting in a system of two ODEs. The purpose of this exercise is to implement the scaled model, test the software, and generalize the model.
a)
Write a function simulate
that can simulate an elastic pendulum using the scaled model.
The function should have the following arguments:
def simulate(
beta=0.9, # dimensionless parameter
Theta=30, # initial angle in degrees
epsilon=0, # initial stretch of wire
num_periods=6, # simulate for num_periods
time_steps_per_period=60, # time step resolution
plot=True, # make plots or not
):
num_periods times \( P \). The time step
is set as \( P \) divided by time_steps_per_period.
The simulate function should return the arrays of
\( x \), \( y \), \( \theta \), and \( t \), where \( \theta = \tan^{-1}(x/(1-y)) \) is
the angular displacement of the elastic pendulum corresponding to the
position \( (x,y) \).
If plot is True, make a plot of \( \bar y(\bar t) \)
versus \( \bar x(\bar t) \), i.e., the physical motion
of the mass at \( (\bar x,\bar y) \). Use the equal aspect ratio on the axis
such that we get a physically correct picture of the motion. Also
make a plot of \( \theta(\bar t) \), where \( \theta \) is measured in degrees.
If \( \Theta < 10 \) degrees, add a plot that compares the solutions of
the scaled, classical, non-elastic pendulum and the elastic pendulum
(\( \theta(t) \)).
Although the mathematics here employs a bar over scaled quantities, the
code should feature plain names x for \( \bar x \), y for \( \bar y \), and
t for \( \bar t \) (rather than x_bar, etc.). These variable names make
the code easier to read and compare with the mathematics.
Hint 1.
Equal aspect ratio is set by plt.gca().set_aspect('equal') in
Matplotlib (import matplotlib.pyplot as plt)
and in SciTools by the command
plt.plot(..., daspect=[1,1,1], daspectmode='equal')
(provided you have done import scitools.std as plt).
Hint 2.
If you want to use Odespy to solve the equations, order the ODEs
like \( \dot \bar x, \bar x, \dot\bar y,\bar y \) such that
odespy.EulerCromer can be applied.
b) Write a test function for testing that \( \Theta=0 \) and \( \epsilon=0 \) gives \( x=y=0 \) for all times.
c) Write another test function for checking that the pure vertical motion of the elastic pendulum is correct. Start with simplifying the ODEs for pure vertical motion and show that \( \bar y(\bar t) \) fulfills a vibration equation with frequency \( \sqrt{\beta/(1-\beta)} \). Set up the exact solution.
Write a test function that
uses this special case to verify the simulate function. There will
be numerical approximation errors present in the results from
simulate so you have to believe in correct results and set a
(low) tolerance that corresponds to the computed maximum error.
Use a small \( \Delta t \) to obtain a small numerical approximation error.
d)
Make a function demo(beta, Theta) for simulating an elastic pendulum with a
given \( \beta \) parameter and initial angle \( \Theta \). Use 600 time steps
per period to get every accurate results, and simulate for 3 periods.
Filename: elastic_pendulum.
This is a continuation Exercise 1.26: Simulate an elastic pendulum.
Air resistance on the body with mass \( m \) can be modeled by the
force \( -\frac{1}{2}\varrho C_D A|\v|\v \),
where \( C_D \) is a drag coefficient (0.2 for a sphere), \( \varrho \)
is the density of air (1.2 \( \hbox{kg }\,{\hbox{m}}^{-3} \)), \( A \) is the
cross section area (\( A=\pi R^2 \) for a sphere, where \( R \) is the radius),
and \( \v \) is the velocity of the body.
Include air resistance in the original model, scale the model,
write a function simulate_drag that is a copy of the simulate
function from Exercise 1.26: Simulate an elastic pendulum, but with the
new ODEs included, and show plots of how air resistance
influences the motion.
Filename: elastic_pendulum_drag.
Test functions are challenging to construct for the problem with
air resistance. You can reuse the tests from
Exercise 1.27: Simulate an elastic pendulum with air resistance for simulate_drag,
but these tests does not verify the new terms arising from air
resistance.
We consider the motion of a planet around a star (the section Two-body gravitational problem). The simplified case where one mass is very much bigger than the other and one object is at rest, results in the scaled ODE model $$ \begin{align*} \ddot x + (x^2 + y^2)^{-3/2}x & = 0,\\ \ddot y + (x^2 + y^2)^{-3/2}y & = 0\tp \end{align*} $$
a) It is easy to show that \( x(t) \) and \( y(t) \) go like sine and cosine functions. Use this idea to derive the exact solution.
b) One believes that a planet may orbit a star for billions of years. We are now interested in how accurate methods we actually need for such calculations. A first task is to determine what the time interval of interest is in scaled units. Take the earth and sun as typical objects and find the characteristic time used in the scaling of the equations (\( t_c = \sqrt{L^3/(mG)} \)), where \( m \) is the mass of the sun, \( L \) is the distance between the sun and the earth, and \( G \) is the gravitational constant. Find the scaled time interval corresponding to one billion years.
c) Solve the equations using 4th-order Runge-Kutta and the Euler-Cromer methods. You may benefit from applying Odespy for this purpose. With each solver, simulate 10000 orbits and print the maximum position error and CPU time as a function of time step. Note that the maximum position error does not necessarily occur at the end of the simulation. The position error achieved with each solver will depend heavily on the size of the time step. Let the time step correspond to 200, 400, 800 and 1600 steps per orbit, respectively. Are the results as expected? Explain briefly. When you develop your program, have in mind that it will be extended with an implementation of the other algorithms (as requested in d) and e) later) and experiments with this algorithm as well.
d) Implement a solver based on the PEFRL method from the section The PEFRL 4th-order accurate algorithm. Verify its 4th-order convergence using an equation \( u'' + u = 0 \).
e) The simulations done previously with the 4th-order Runge-Kutta and Euler-Cromer are now to be repeated with the PEFRL solver, so the code must be extended accordingly. Then run the simulations and comment on the performance of PEFRL compared to the other two.
f) Use the PEFRL solver to simulate 100000 orbits with a fixed time step corresponding to 1600 steps per period. Record the maximum error within each subsequent group of 1000 orbits. Plot these errors and fit (least squares) a mathematical function to the data. Print also the total CPU time spent for all 100000 orbits.
Now, predict the error and required CPU time for a simulation of 1 billion years (orbits). Is it feasible on today's computers to simulate the planetary motion for one billion years?
Filename: vib_PEFRL.
This exercise investigates whether it is feasible to predict planetary motion for the life time of a solar system.