Page 1 of 9
Assignment 9: Explicit Time-Marching Schemes
Due Date: Wednesday 29 April 11:59pm.
Overview
This assignment brings together object-oriented programing with numerical time-marching
methods. For the first part of the assignment you will implement some explicit time-marching
schemes within an object-oriented framework to solve systems of ODEs that model initial value
problems. For the second part of the assignment you will model the dynamics of a crane cable
that might fail when the crane mechanism seizes up while lowering a load.
Instructions and Deliverables
Part 1: ode module
Create a module file named ode.py. This file must define the following classes:
- ODESolver: a base class for all time-marching methods.
- ForwardEuler: a subclass of ODESolver that implements forward Euler.
- RK4: a subclass of ODESolver that implements the 4th-order Runge-Kutta method.
Before describing the methods to implement for each of the above classes, we first illustrate how
the classes would be used to solve an initial-value problem (IVP). Suppose we want to solve the
bacteria concentration IVP described in the lectures:
This could be solved and the solution displayed using the ForwardEuler class with the
following code:
import ode
import numpy as np
def bacteria(c, t):
"""Returns the ODE right-hand-side function
for the bacteria decontamination model"""
f = -0.6*c
return f
tspan = np.array([0.0, 7.0])
num_steps = 10
y0 = np.array([1.0e7]) # note that the IC is passed as an array
solver = ode.ForwardEuler(bacteria) # construct variable of type ForwardEuler
t, y = solver.solve_ivp(tspan, y0, num_steps) # uses method solve_ivp
Page 2 of 9
print()
print(" i time count ")
print("--- ---- --------")
for i in range(y.shape[0]):
string = '{:3d}'.format(i ) + " "
string += '{:4.2f}'.format(t[i] ) + " "
string += '{:8.2e}'.format(y[0, i])
print(string)
Look closely at the follwing two lines in particular. You must make yourself familiar with this
object oriented approach to applying the numerical methods to understand this assignment.
solver = ode.ForwardEuler(bacteria) # construct variable of type ForwardEuler
– ode is the imported module file that contains various numerical methods for solving
ODEs.
– solver is a variable constructed to be of type ForwardEuler, a class that is
defined in module ode. solver gets its properties from class ForwardEuler
(and its super class ODESolver).
– bacteria is a function (defined above) attached to solver via the __init__
function that ForwardEuler will inherit from super class ODESolver.
Through this construction, variable solver now carries bacteria with it to use in the
solution method.
t, c = solver.solve_ivp(tspan, c0, num_steps) # uses method solve_ivp
– solver is the variable constructed above with the function bacteria attached.
– solve_ivp is a method in class ForwardEuler (that it actually inherits from
super class ODESolver) which will apply the forward Euler numerical method
to solve the initial value problem. solve_ivp takes three arguments:
0
B tspan is a two element NumPy array holding tspan[0] = t and
tspan[1] f
= t .
B y0 is a 1D NumPy array with as many elements as there are ODEs in
the system of ODEs to be solved. As shown above, the bacteria
example has only one ODE so a one-element array holds the one
initial condtion.
0 f B num_steps is the number of time steps to use between t and t .
To solve the same IVP with the 4th-order Runge-Kutta method, we would only need to replace
the ForwardEuler constructor with the RK4 constructor.
Page 3 of 9
The ODESolver base class shall meet the following requirements:
C The class constructor will take in a function and define the _func member based on
this function. Thus, the first-line definition should be
def __init__(self, func):
For any variable constructed by ODESolver we want the argument passed to
__init__ to be “attached” to that constructed variable.
Recall from Lecture 16 how __init__ works. In the example on pp. 6-11,
when we construct a variable of type Vector we associate parameters in
function __init__ (x and y) by which arguments can be assigned to the
constructed variable (self._x and self._y). In ODESolver we do the
same, except that we associate a parameter that can be a function: func, in the
same way that we can pass functions to other functions (e.g., see Lecture 19 p. 14
where the forward_euler function has parameter func). Note that we want
to use an underscore as the first character to signal that the attached function is
private: that when interacting with your class no one should arbitrarily reassign
self._func for any reason. (The same way one should not write something
like math.pi = 3; perhaps it should instead be named math._pi to signal
that.)
C The class should include the member function solve_ivp with the following
first-line definition:
def solve_ivp(self, tspan, y0, num_steps):
This was described in detail above.
As with the in-class exercises, raise a ValueError with a useful message if the
arguments to tspan or num_steps contain an obvious error.
The method solve_ivp should return two arrays, t and y, in that order. The
array t should hold all num_steps+1 time values including tspan[0] and
tspan[1]. The array y should be a two-dimensional NumPy array with
y.shape[0] = y0.size and y.shape[1] = t.size; that is, the
number of “rows” in y should equal the number of dependent variables (and, thus,
the number of ODEs in the system), and the number of “columns” in y should
equal the number of values in t.
Finally, solve_ivp should not be re-implemented in the subclasses. The idea
is to write this routine once and use it for all the subclasses via inheritance. The
Page 4 of 9
solution will thus be calculated by solve_ivp as
for n in range(0, num_steps):
y[:,n+1] = self.step(t[n], dt, y[:,n])
Compare this to how solutions were calculated in previous in-class exercises (e.g.,
using forward Euler, and note two changes:
– the solution variable, y, is now a 2D NumPy array
– the step method comes from self. For the bacterial colony example
above, self is the forwardEuler type that solver was
constructed to be. Thus, self.step will be the same as calling
forwardEuler.step – the step function defined in the
forwardEuler sub-class. (If the user instead chose to construct
solver to be of type RK4, then this would be the same as calling
RK4.step and would invoke the step function defined in that
sub-class.)
The ODESolver class is an abstract base class: it should be used to derive new
classes such as ForwardEuler and RK4, and it should not be used to solve
IVPs on its own. To enforce this intent, provide ODESolver with a step method
with the following first-line definition:
def step(self, t_old, dt, y_old):
However, for the ODESolver base class only, this function should only raise a
NotImplementedError and not do anything else.
The ForwardEuler subclass should meet the following requirements.
C The ForwardEuler class should be derived from ODESolver.
C The ForwardEuler class should implement only one method, namely step.
The step method should have the identical first-line definition as for super class
ODESolver:
def step(self, t_old, dt, y_old):
where t_old n
is the time value at the current time (referred to as t in the lecture
slides), dt is the time step size )t, and y_old is the variable value at the current
n n+1 time (y ). step will return the value y as defined by the forward Euler
numerical method.
Page 5 of 9
The RK4 subclass should meet the following requirements.
C The RK4 class should be derived from ODESolver.
C The RK4 class should implement only one method, namely step. The step method
should have the first-line definition shown above for forwardEuler, and it
n+1 should return the value y as defined by the classical 4th-order Runge-Kutta
method.
As usual, provide docstrings for the ode module, the classes, and all methods.
Part 2: cable mechanism seizure failure
Background and Models
Consider the classic crane from IEA, pictured at right. A force
balance on the load when static gives:
where
+y is in the downward (lowering) direction
W is the weight of the load [N]
T is the tension in the cable [N]
m is the mass of the load [kg] and g is the constant of gravitation [9.80665 m/s ]
2
Recall that “static” means non-accelerating, so when a load is being lowered at a constant
velocity it is also in static equilibrium. From this model, a cable that can support a 7 kN tension
can theoretically lower a 7 kN load at constant velocity.
However, what if the cable mechanism suddenly seizes while lowering the load? Suppose that
the load is being lowered at some velocity v and suddenly the cable mechanism stops. From
Physics I, the dynamics of the load involves:
– inertia (the “ma” in F = ma),
– a damping/restoring force (not present in this problem as described so far) proportional
to the dynamic velocity of the load,
– the elastic tension in the cable (a spring force) due to elongation from the change in
momentum of the load when the cable seizes,
– the static tension in the cable, and
– the weight of the load,
as shown in this equation:
Page 6 of 9
where
a(t) is the acceleration of the load [m/s ]
2
v(t) is the velocity of the load [m/s]
y(t) is the displacement of the load from the point it stops [m]
b is a damping coefficient [N@s/m]
k is the spring constant of the elastic cable (stretching from the inertia) [N/m]
The initial conditions are that the initial position y(0) is arbitrarily set to zero, and the initial
velocity v(0) is the ve initial locity at which the load is being lowered, v .
From the static model W ! T = 0. Noting that dy(t)/dt = v(t) and d y(t)/dt = a(t), and 2 2
multiplying through by !1, the equation of motion of the dynamic load can be written as:
After some algebra this can also be written as a system of first order ODEs:
In Physics I and possibly in Introduction to Differential Equations you have already seen these
derived analytically, which is useful for checking a numerical solution. For a review of the
theory of harmonic oscillators you may find the following helpful:
https://physics.info/sho/
http://hyperphysics.phy-astr.gsu.edu/hbase/shm.html#c1
https://en.wikipedia.org/wiki/Harmonic_oscillator
However, you may not use the analytical solutions for this assignment (except as optional
checks). All results should be computed from your numerical solution.
Frequency f of undamped (b = 0) harmonic oscillator
Note: if you plot y(t) vs t for 1 second you should be able to count f periods if
you’ve solved the system correctly.
Page 7 of 9
Max max imum tension T in the cable
Max max imum stress F in the cable
Criterion for cable failure
Deliverables
In a module file titled hw9.py include the following (as well as the usual useful doctstrings and
comments throughout).
final C all model input data (e.g., mass, spring constant, etc.) and numerics input (t , number
of steps) should be included as module variables near the top of your file that you
can easily adjust. Include the units for the variables as in-line comments. For
example,
g = 9.80665 # constant of gravitation [m/s**2]
You may find it helpful to optionally print this out for your information; e.g.,
print()
print("LOAD ----- ----- -----")
print(" mass, m [kg] = ", '{:7.1f}'.format(m))
print(" weight, W [N] = ", '{:7.1f}'.format(W))
print("CABLE ---- ----- -----")
print(" spring constant, k [N/m] = ", '{:12.2e}'.format(k))
# etc.
C a function that defines the right hand side of the ODE model to solve.
This function will input the array of state variables and the time and output the
right hand side of the system of ODEs model; e.g.,
Page 8 of 9
def crane_cable(state, t):
'''Returns the ODE right-hand-side function
for the crane cable mechanism seizure failure model
Input: state[:,i] - 1D slice of state variables at i-th time step
t - time at i-th time step
Returns: dsdt[:,i] - 1D ODE rhs function evaluated at input state
'''
C code to import and apply an appropriate ODE solver from your ode.py module.
C code to create a professional .png plot of load position versus time for 1 second. As
described above, use this plot to verify that you have solved the model correctly.
Use the if(Make_Plot): technique described in Assignment 8 to hide any
references to matplotlib from Submitty before you submit file hw9.py.
C code to compute the principle “underlined” quantities listed above and print them out
in a well formatted table that includes units. When you are ready to submit your
assignment to Submitty, copy this output into comments at the end of file hw9.py.
There will be two sets of output, described subsequently. The following is an
example of how you might do this. Please output the values in the units shown.
print("RESULTS ** ***** ***** ***** *****")
print(" frequency =", '{:6.2f}'.format(frequency) , "[Hz]")
print(" amplitude =", '{:6.2f}'.format(max_amplitude*1.0e3), "[mm]")
print(" tension =", '{:6.2f}'.format(max_tension/1.0e3) , "[kN]")
print(" stress =", '{:6.2f}'.format(max_stress /1.0e6) , "[MPa]")
C a comment block (e.g., text with000 A)... B)... C)... D)...000
at the very end of your code briefly answering the following questions.
A) How do you know you’ve used an appropriate ODE solver?
[Hint: look over your notes that describe the stability of the methods you
have available in your ode.py module.]
B) State a criterion for deciding how many time steps to use (i.e., how to get a
small enough )t) and state the number of time steps you used to meet that
criterion. You do not necessarily need a quantitative criterion.
[Hint: our derivations in the notes were about the stability of the solution,
not the accuracy. How can you judge accuracy? Consider the
convergence criteria we’ve studied previously this semester, and run
Page 9 of 9
numerical experiments with various numbers of steps to come to a
conclusion.]
C) Choose a unique cable strength between 50.0 and 100 MPa. Use your code to
manually (or automatically, if you like) compute the diameter of the cable
(in mm, to two significant figures rounded up) needed to survive this
accident with a factor of safety of 3.0. (Note that factors of safety are still
commonly used in design scoping work, but are now deprecated in many
safety related engineering analyses.)
[Hint: “manually” means you can manually adjust an input value and
rerun your code to get a desired output. Ideally you would implement a
root-finding procedure to do this, but unfortunately (?) we believe the time
remaining in this semester is too short to ask for that.]
D) A vendor has a damping bar attachment that can mitigate a cable seizure
accident with a damping coefficient of b = 8000.0 N@s/m. For the cable
strength you chose in Question C, qualitatively how does this change the
cable diameter – a lot or a little, and why?
[Hint: it is recommended that you carefully examine and compare the two
plots and recall how the maximum stress is calculated to consider this
question.]
Please note that you have a lot of freedom in how you code and answer this assignment. The
point is not to get a “right” answer, it is to identify a safe design that you are confident it without
knowing exactly what the “right” answer is.
Use the following data in your model. Recall that these should apper as module variables.
mass of load [kg] = 500.0 kg
damping coefficient [N@s/m] = 0.0 N.s/m (except for question D)
spring constant of cable in N/m = 1.05@10 N/m
6
constant of gravitation = 9.80665 m/s
cable speed = 0.1 m/s (This assignment should be
informative about why the speed is
generally so slow; experiment with
other values to see the effect.)
final t = 1.0 s