1

The question body:

A skydiver of mass m in a vertical free fall experiences an aerodynamic drag force F=cy'² ('c times y prime square') where y is measured downward from the start of the fall, and y is a function of time (y' denotes the derivative of y w.r.t time). The differential equation describing the fall is:

y''=g-(c/m)y'²

(where g = 9.80665 m/s^2; c = 0.2028 kg/m; m = 80 kg). And y(0)=y'(0)=0 as this is a free fall.

Task: The function must return the time of a fall of x meters, where x is the parameter of the function. The values of g, c and m are given below.

The Runge-Kutta function is defined as follows:

from numpy import *
def runge_kutta_4(F, x0, y0, x, h):
   '''
   Return y(x) given the following initial value problem:
   y' = F(x, y)
   y(x0) = y0 # initial conditions
   h is the increment of x used in integration
   F = [y'[0], y'[1], ..., y'[n-1]]
   y = [y[0], y[1], ..., y[n-1]]
   '''
   X = []
   Y = []
   X.append(x0)
   Y.append(y0)
   while x0 < x:
       k0 = F(x0, y0)
       k1 = F(x0 + h / 2.0, y0 + h / 2.0 * k0)
       k2 = F(x0 + h / 2.0, y0 + h / 2 * k1)
       k3 = F(x0 + h, y0 + h * k2)
       y0 = y0 + h / 6.0 * (k0 + 2 * k1 + 2.0 * k2 + k3)
       x0 += h
       X.append(x0)
       Y.append(y0)
   return array(X), array(Y)

And this is what I've done so far:

def prob_1_8(x)
    g = 9.80665  # m/s**2
    c = 0.2028  # kg/m
    m = 80  # kg

    def F(x, y):
        return array([
            y[1],
            g - (c / m) * ((y[1]) ** 2)
        ])

    X, Y = runge_kutta_4(F, 0, array([0, 0]), 5000, 1000)
    for i in range(len(X)):
        if X[i] == 5000:
            return Y[i]

However, when I tried to print prob_1_8(5000), the number looks ridiculous and it displayed:

RuntimeWarning: overflow encountered in double_scalars.

According to the answer provided, I should get a value close to 84.8 when x=5000. Can someone help me with this? I don't know what's the problem and how to fix it.

Community
  • 1
  • 1
Ziming Wang
  • 25
  • 3

1 Answers1

0

Please contemplate the function call of X, Y = runge_kutta_4(F, 0, array([0, 0]), 5000, 1000). You are integrating over a time span of 5000 sec > 1 hour in steps of 1000 sec > 16 min. It is intuitively clear that this will be imprecise, as most of the acceleration will happen in the first 10 sec.

Then the question is what exactly you are trying to filter out with the loop. Is it the speed after this time?

The limit speed is where the right side is zero, at vmax=sqrt(g*m/c) = 62.1972 = 223.91 km/h, the claimed value of 84.8 can not be reached as a speed starting from rest. The fall time to a distance of x will be a little more than x/vmax, so you could use tmax = 100+x/vmax in T, Y = runge_kutta_4(F, t0, y0, tmax, 1).

Integrating in 1 sec time steps and looking for the speed after 5000 meters fall distance gives a result of 85 sec, distance 5013.33465614 m, speed 62.1972 m/s which is as to be expected close to the limit speed.

You can get a more precise time value by using (reverse) linear interpolation, then at time about 84.786 sec you reach distance 5000 m with a speed 62.1972 m/s. This again is compatible with the claimed result value, which now is a time, not a velocity.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51