4. Numerics

A marble is thrown into a transparent container filled with a viscous liquid.  Snapshots of its trajectory are taken at regular time intervals to record a set of samples $$z_i = z(i\Delta t)$$ of its vertical position. Δt is the time interval between the samples.

Under the assumption that the friction force is linearly dependent on the marble's velocity, the (vertical) velocity and position should have the following form: 

$$v_t = (v_0 + \frac{g}{\alpha}) e^{-\alpha t} - \frac{g}{\alpha}$$

$$z_t = z_0 + \frac{1}{\alpha} (v_0+\frac{g}{\alpha})(1-e^{-\alpha t}) - \frac{g}{\alpha} t$$

where

$$ z_0 \mbox{ is the initial vertical position,} $$

$$ v_0 \mbox{ is the initial velocity in the vertical direction,} $$

$$ g \mbox{ is the gravitation accelleration, and} $$

and

$$ \alpha \mbox{ is the friction rate.} $$

The friction rate α is the propery we are after here. The time between samples, Δt, is known but all other parameters should also be considered to be unknown.  (Note: you do not have to take any horizontal motion into account.)

Knowing these forms make it possible to estimate the friction force from the samples using a module called friction.  It contains two functions:

 double frictionrate(double dt, const rvector<double>& v);
 rvector<double> numdiff(double dt, const rvector<double>& z);

The first function, frictionrate, computes the friction rate based on velocity samples v taken a time dt apart.

The second function, numdiff, is needed because the first function requires velocity samples, which the experiment does not provide.  We can, however, estimate the velocities using finite differences of the position samples.  This is what the numdiff function does. 

So given a set of samples in an rvector z, the friction rate can be estimated with "rate = frictionrate(dt,numdiff(dt,z));".

For the assignment you should

• Write at least 2 unit tests for the numdiff function and two for the frictionrate function.  These tests should create test data for which the result is known, and then test their output.

    For at least one of the tests for each of these two functions, use the equations above to generate this test data, with the following parameters: α=0.125, v₀=10, z₀=0, g=9.8, Δt=0.25, and t ranging from 0 and 100.

    For the unit tests, use one of the following frameworks: catch2, Boost.Test, or googletest. All three are available as modules on the Teach cluster ("module load catch boost googletest"). Catch2 is perhaps the simplest to get started with, but the lecture slides on software testing have examples for Boost.Test.

• Find and build a test case for which the frictionrate function fails. It may help to know that the frictionrate function actually computes the following:

$$\alpha = \frac{1}{(n-2)\Delta t}\sum_{i=1}^{n-2} \ln \frac{v_{i+2}-v_{i+1}}{v_{i+1}-v_{i}}$$

Correction (March 1): That is not what is actually calculated in the code, instead, the function frictionrate computes the following: 

$$\alpha = \frac{-1}{\Delta t}\ln \left[\frac{1}{n-2}\sum_{i=1}^{n-2} \frac{v_{i+2}-v_{i+1}}{v_{i+1}-v_{i}} \right]$$

and this isn't particularly robust and can get ``nan''s or ``inf''s as output fairly easily. 

• Can you think of a way to modify the frictionrate function to be more robust?
• Write a separate program to assess the discretization error as  Δt increases, i.e., it should generate test data for increasing Δt for the range of t, and use it in the determination of α.  It should show that small Δt's give accurate results but large Δt's do not. 

Use a makefile and git throughout your assignment, i.e. whenever you write a test, create a target for it in the Makefile, such that "make TARGETTEST" runs the test called TARGETTEST and commit. Also create an overall "test" target in your Makefile that runs all of tests. We expect to see several meaningful git commits.  Make sure to commit all necessary files to run the code and its tests. Upload all code, the git log, and a text file with a short answer to the robustness question above.