GEOL 5470    Humphrey 2022

Homework, for Feb 8 2022

I am assuming you have got your code for the falling rock problem running, hopefully both the SciPy ODEINT version and also your own version using a simple Euler step technique. You have run your code using lots of small steps, what we want to see is if taking big steps introduces errors. So run your falling rock code with only 2 or 3 time steps from 0 all the way to terminal velocity.  So we are consistent, use an aggregate ‘c’ value ( in g – c * v² terms ) of about .02, that will mean your code must span about 7 to 10 secs.

1. Run both your Euler code and the Scipy odeint code.  Can you explain why your Euler code doesn’t work very well for long time steps?

2. Write a Runge-Kutta routine, as discussed in class, and as shown in the note below, instead of your simple “Euler” code for the prediction of each step.  Does this help? 

3. Plot the results of a. your simple Euler code, b. scipy odeint code, and finally c. for a Runge Kutta  version.  Experiment with your programs with different time steps. Which program handles the big time steps the best?

A note on the Runge-Kutta method.  I want you to code up the Runge-Kutta (this is the almost universal 4^th order method). The main difficulty is with keeping track of the intermediate values. Note you have to calculate the intermediate values in the order shown.  Here is the scheme to calculate one time step for the velocity, if we call the derivative ( [ g – c * v² ] ) with the symbol f.

v_i+1 = v_i + delt * [ 1/6  f(v_i, t_i) + 1/3  f(v^*_i+1/2 , t_i+1/2) + 1/3   f(v^**_i+1/2, t_i+1/2) + 1/6   f(v^*_i+1, t_i+1) ]

where these are the intermediate steps (in order):

1: v^*_i+1/2 = v_i + delt/2   f(v_i , t_i)

2: v^**_i+1/2 = v_i + delt/2   f(v^*_i+1/2  , t_i)

3: v^*_i+1   =   v_i + delt/2   f(v^**_i+1/2  , t_i)

Note you actually need to calculate all the starred ‘v’ values first (in order) to plug into the final equation for the prediction of the final v_{i+1 .}  And notice that instead of just one call to your function for ‘f’, there are 6 calls for one time step.

It is best to write the R-K process as a function that you call with the value of v_i and t, and which will return the value of v_{i+1 .}  The R-K function will call your ‘f’ function (your derivative function).