python - 如何使用 python 中提供的 Runge-Kutta 方法解决以下问题

Question

问题正文:

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.

龙格-库塔函数定义如下:

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)

这就是我到目前为止所做的:

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]

但是，当我尝试打印 prob_1_8(5000) 时，这个数字看起来很荒谬，并且显示:

RuntimeWarning: overflow encountered in double_scalars. 

根据提供的答案，当 x=5000 时，我应该得到接近 84.8 的值。有人可以帮我弄这个吗？我不知道问题是什么以及如何解决。

Answer 1

请考虑X, Y = runge_kutta_4(F, 0, array([0, 0]), 5000, 1000)的函数调用。您正在积分的时间跨度为 5000 sec > 1 hour步长为1000 sec > 16 min 。直观上很明显，这将是不精确的，因为大部分加速将发生在前 10 秒内。

<小时/>

那么问题是你到底想用循环过滤掉什么。是这次之后的速度吗？

<小时/>

极限速度是右侧为零的位置，即 vmax=sqrt(g*m/c) = 62.1972 = 223.91 km/h ，声称的值84.8无法达到从静止开始的速度。下降到距离x的时间会比 x/vmax 多一点，所以你可以使用 tmax = 100+x/vmax在T, Y = runge_kutta_4(F, t0, y0, tmax, 1) .

<小时/>

以 1 秒时间步长积分并查找 5000 米坠落距离后的速度，得出的结果为 85 sec ，距离5013.33465614 m ，速度62.1972 m/s这是预期的接近极限速度。

您可以通过使用(反向)线性插值获得更精确的时间值，然后在大约 84.786 sec 的时间您到达距离5000 m以速度62.1972 m/s 。这再次与声称的结果值兼容，现在是时间，而不是速度。