Тело вопроса:
Парашютист массой м при вертикальном свободном падении испытывает силу аэродинамического сопротивления F = cy'² ('c раз y < / em> prime square '), где y измеряется вниз от начала падения, а y является функцией времени (y' обозначает производную от y по времени). Дифференциальное уравнение, описывающее падение:
y''=g-(c/m)y'²
(где g = 9,80665 м / с ^ 2; c = 0,2028 кг / м; m = 80 кг). И y (0) = y '(0) = 0, так как это свободное падение.
Задача: функция должна возвращать время падения x метров, где x - параметр функции. Значения g, c и m приведены ниже.
Функция Рунге-Кутта определяется следующим образом:
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
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)
RuntimeWarning: overflow encountered in double_scalars.
Согласно предоставленному ответу, я должен получить значение, близкое к 84,8, при x = 5000. Может кто-то помочь мне с этим? Не знаю, в чем проблема и как ее исправить.
