Реализация начальных условий для численно решаемого дифференциального уравнения

Представьте, что кто-то прыгает с балкона под определенным углом. theta и скорость v0 (высота балкона обозначается как ystar). Глядя на эту проблему в 2D и рассматривая перетаскивание, вы получаете систему дифференциальных уравнений, которая может быть решена с помощью метода Рунге-Кутты (я выбираю явную среднюю точку, не уверенный, какова таблица мясника для этого). Я реализовал это, и он отлично работает, для некоторых начальных условий я получаю траекторию движущейся частицы.

Моя проблема в том, что я хочу исправить два начальных условия (начальная точка на оси x равна нулю, а на оси y - ystar) и убедитесь, что траектория проходит через определенную точку на оси x (назовем ее xstar). Для этого, конечно, существует множество комбинаций двух других начальных условий, которые в этом случае являются скоростями в направлении x и y. Проблема в том, что я не знаю, как это реализовать.

Код, который я использовал для решения проблемы до этого момента:

1) Реализация метода Рунге-Кутты

import numpy as np 
import matplotlib.pyplot as plt

def integrate(methode_step, rhs, y0, T, N):
    star = (int(N+1),y0.size)
    y= np.empty(star)
    t0, dt = 0, 1.* T/N
    y[0,...] = y0
    for i in range(0,int(N)):
        y[i+1,...]=methode_step(rhs,y[i,...], t0+i*dt, dt)
    t = np.arange(N+1) * dt
    return t,y

def explicit_midpoint_step(rhs, y0, t0, dt):
    return y0 + dt * rhs(t0+0.5*dt,y0+0.5*dt*rhs(t0,y0))

def explicit_midpoint(rhs,y0,T,N):
    return integrate(explicit_midpoint_step,rhs,y0,T,N)

2) Реализация правой части дифференциального уравнения и вспомогательных параметров

A = 1.9/2.
cw = 0.78
rho = 1.293
g = 9.81

# Mass and referece length
l = 1.95
m = 118

# Position
xstar = 8*l
ystar = 4*l

def rhs(t,y):
    lam = cw * A * rho /(2 * m)
    return np.array([y[1],-lam*y[1]*np.sqrt(y[1]**2+y[3]**2),y[3],-lam*y[3]*np.sqrt(y[1]**2+y[3]**2)-g])

3) решение проблемы с ним

# Parametrize the two dimensional velocity with an angle theta and speed v0

v0 = 30
theta = np.pi/6

v0x = v0 * np.cos(theta)
v0y = v0 * np.sin(theta)

# Initial condintions
z0 = np.array([0, v0x, ystar, v0y])

# Calculate solution
t, z = explicit_midpoint(rhs, z0, 5, 1000)

4) Визуализация

plt.figure()
plt.plot(0,ystar,"ro")
plt.plot(x,0,"ro")
plt.plot(z[:,0],z[:,1])
plt.grid(True)
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.show()

Чтобы конкретизировать вопрос: с учетом этого, как мне найти все возможные комбинации v0 и theta так, чтобы z[some_element,0]==xstar

Я, конечно, пробовал кое-что, в основном метод грубой силы, чтобы исправить theta и затем опробовать все возможные скорости (в это имеет смысл) но, наконец, не знал, как сравнить полученные массивы с желаемым результатом ...

Поскольку это, в основном, проблема с кодированием, я надеюсь, что переполнение стека - подходящее место для обращения за помощью ...

EDIT: В соответствии с просьбой, вот моя попытка решить проблему (замена 3) и 4) сверху).

theta = np.pi/4.
xy = np.zeros((50,1001,2))
z1 = np.zeros((1001,2))
count=0

for v0 in range(0,50):
    v0x = v0 * np.cos(theta)
    v0y = v0 * np.sin(theta)
    z0 = np.array([0, v0x, ystar, v0y])

    # Calculate solution
    t, z = explicit_midpoint(rhs, z0, 5, 1000)    

    if np.around(z[:,0],3).any() == round(xstar,3):
        z1[:,0] = z[:,0]
        z1[:,1] = z[:,2]
        break
    else:
        xy[count,:,0] = z[:,0]
        xy[count,:,1] = z[:,2]
        count+=1


plt.figure()
plt.plot(0,ystar,"ro")
plt.plot(xstar,0,"ro")

for k in range(0,50):    
    plt.plot(xy[k,:,0],xy[k,:,1])

plt.plot(z[:,0],z[:,1])
plt.grid(True)
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.show()

Я уверен, что я неправильно использую функцию .any(), идея заключается в том, чтобы округлить значения z[:,0] до трех цифр и сравнить их с xstar, если он совпадает, цикл должен завершиться и повторно запустить текущий z, если нет, он должен сохранить его в другом массиве, а затем увеличить v0.

4 голоса | спросил Sito 14 J000000Saturday18 2018, 01:01:30

3 ответа


0
Изменить 2018-07-16Здесь я выкладываю исправленный ответ с учетом сопротивления воздуха.Ниже приведен скрипт на python для вычисления набора значений ---- +: = 0 =: + ---- так, чтобы воздушная траектория проходила через ---- +: = 1 =: + ----в какое-то время ---- +: = 2 =: + ---- .Я использовал траекторию без аэродинамического сопротивления в качестве первоначальной догадки, а также чтобы угадать зависимость ---- +: = 3 =: + ---- от ---- +: = 4 =: + ----для первой доработки.Количество итераций, необходимых для получения правильного значения ---- +: = 5 =: + ----, обычно составляло от 3 до 4. Сценарий на моем ноутбуке завершился за 0,99 секунды, включая время генерации фигур.Скрипт генерирует две фигуры и один текстовый файл.Черные точки указывают набор решений ---- +: = 7 =: + ----Желтая линия указывает на ссылку ---- +: = 8 =: + ----, которая была бы решением, если бы не было воздушного сопротивления.Проверка того, что траектория (синяя сплошная линия) проходит через ---- +: = 10 =: + ----, когда ---- +: = 11 =: + ---- , выбрана из набора решений.Черная пунктирная линия показывает траекторию без сопротивления воздухом в качестве ориентирасодержит числовые данные ---- +: = 13 =: + ----, а также время посадки ---- +: = 14 =: + ---- и количество итераций, необходимых для поиска -- +: = 15 =: + ---- .Здесь начинается сценарий.Вот рисунок ---- +: = 17 =: + ---- .fig_xdrop_v0_theta.pngВот рисунок ---- +: = 18 =: + ---- .fig_traj_sample.pngИзменить 2018-07-15Я понял, что упустил из виду, что вопрос касается сопротивления воздуха.Какой позор для меня.Итак, мой ответ ниже не является правильным.Я боюсь, что удаление моего ответа само по себе похоже на сокрытие ошибки, и я оставляю это ниже на данный момент.Если люди думают, что раздражает неправильный ответ, я в порядке, кто-то удалит его.Дифференциальное уравнение на самом деле может быть решено вручную, и оно не требует больших вычислительных ресурсов, чтобы определить, как далеко человек достигает балкона от земли в зависимости от начальной скорости ---- +: = 19 =: +---- и угол ---- +: = 20 =: + ---- .Затем вы можете выбрать условие ---- +: = 21 =: + ---- так , чтобы ---- +: = 22 =: + ---- из этой таблицы данных.Запишем горизонтальные и вертикальные координаты человека во время ---- +: = 23 =: + ---- is ---- +: = 24 =: + ---- и ---- +:= 25 =: + ---- соответственно.Я думаю, что вы взяли ---- +: = 26 =: + ---- у стены здания и ---- +: = 27 =: + ---- за уровень земли, и я так и делаюздесь тоже.Скажем, горизонтальная и вертикальная скорость человека во время ---- +: = 28 =: + ---- являются ---- +: = 29 =: + ---- и ---- +:= 30 =: + ---- соответственно.Начальные условия в ---- +: = 31 =: + ---- задаются какУравнение Ньютона, которое вы решаете,где ---- +: = 34 =: + ---- - масса человека, а ---- +: = 35 =: + ---- - константа, которую я не знаю по-английскиимя, но мы все знаем, что это такое.Из уравнения(3),Используя это с уравнением.(1),где мы также использовали начальныйусловие.---- +: = 38 =: + ---- означает целое число от ---- +: = 39 =: + ---- до ---- +: = 40 =: + ---- .Из уравнения(4),где мы использовали начальное условие.Используя это с уравнением.(3) а также с использованием начального условия,где ---- +: = 43 =: + ---- означает ---- +: = 44 =: + ---- квадрат.Из выражения для ---- +: = 45 =: + ---- мы можем получить время ---- +: = 46 =: + ----, когда человек падает на землю.То есть ---- +: = 47 =: + ---- .Это уравнение может быть решено квадратичной формулой (или чем-то похожим именем) какгде я использовал условие ---- +: = 49 =: + ---- .Теперь мы знаем расстояние от балкона, которого достиг человек, когда он упал на землю, как ---- +: = 50 =: + ---- .Используя выражение для ---- +: = 51 =: + ---- выше,На самом деле ---- +: = 53 =: + ---- зависит от ---- +: = 54 =: + ---- и ---- +: = 55 =: + ---- кака также ---- +: = 56 =: + ---- и ---- +: = 57 =: + ---- .Вы держите в качестве констант ---- +: = 58 =: + ---- и ---- +: = 59 =: + ---- и хотите найти все ---- +: = 60=: + ---- такое, что ---- +: = 61 =: + ---- для заданного ---- +: = 62 =: + ---- значение.С правой стороны уравнения.(5) можно вычислить дешево, вы можете настроить сетки для ---- +: = 63 =: + ---- и ---- +: = 64 =: + ---- и вычислить ---- +: = 65 =: + ---- на этой 2D сетке.Затем вы можете увидеть, где примерно находится множество решений ---- +: = 66 =: + ---- лежит.Если вам нужно точное решение, вы можете выбрать сегмент, который включает решение из этой таблицы данных.Ниже приведен скрипт на python, демонстрирующий эту идею.Я также прикрепляю сюда рисунок, сгенерированный этим скриптом.Желтая кривая - это решение, установленное ---- +: = 67 =: + ---- таким образом, что человек упал на землю в ---- +: = 68 =: + ---- от стены, когда ---- +: = 69 =: + ---- и ---- +: = 70 =: + ----, как вы установили.Синяя цветовая координата указывает ---- +: = 71 =: + ---- , т. Е. Как далеко человек прыгнул с балкона по горизонтали.Обратите внимание, что в данном ---- +: = 72 =: + ---- (выше порогового значения aruond ---- +: = 73 =: + ---- ), есть два -- +: = 74 =: + ---- значения (два направления для того, чтобы человек проецировал себя) для достижения намеченной точки ---- +: = 75 =: + ---- .Меньшая ветвь значения ---- +: = 76 =: + ---- может быть отрицательной, что означает, что человек может прыгнуть вниз, чтобы достичь намеченной точки, пока начальная скорость достаточно высока.Сценарий также генерирует файл данных ---- +: = 77 =: + ---- , который содержит сегменты, включающие решение.fig_xdrop_v0_theta.png
ответил norio 15 J000000Sunday18 2018, 09:10:08
0
Поэтому после некоторой попытки я нашел способ добиться того, чего хотел ... Это метод грубой силы, о котором я упоминал в своем начальном посте, но, по крайней мере, теперь он работает ...Идея довольно проста: определить функцию ---- +: = 0 =: + ----, которая находит для данного ---- +: = 1 =: + ---- a ---- +: = 2 =: + ---- .В этой функции вы берете начальное значение для ---- +: = 3 =: + ---- (я выбираю 8, но это было только предположение от меня), затем берете начальное значение и проверяете с помощью ---- +: = 4 =: + ---- определяет, насколько далеко интересная точка находится от ---- +: = 5 =: + ---- .Интересная точка в этом случае может быть определена путем установки всех точек на оси x, которые больше чем ---- +: = 6 =: + ---- в ноль (и их соответствующие значения y), а затем обрезкаиз всех нулей с ---- +: = 7 =: + ---- , теперь последний элемент соответствует желаемому результату.Если выходной сигнал разностной функции меньше критического значения (в моем случае 0,1), передайте ток ---- +: = 8 =: + ---- , если нет, увеличьте его на ---- +: = 9 =: + ---- и сделать то же самое снова.Код для этого выглядит следующим образом (снова заменяя 3) и 4)):Проблема с этим состоит в том, что для вычисления требуется вечность (с текущими настройками около 5-6 минут).Если у кого-то есть подсказки, как улучшить код, чтобы он стал немного быстрее, или есть другой подход, он все равно будет оценен.
ответил Sito 14 J000000Saturday18 2018, 16:54:37
0
Предполагая, что скорость в направлении x никогда не падает до нуля, вы можете принять x в качестве независимого параметра вместо времени.Тогда вектор состояния - это время, положение, скорость, и векторное поле в этом пространстве состояний масштабируется так, чтобы компонент vx всегда был равен 1. Затем интегрируем от нуля до xstar, чтобы вычислить состояние (приближение), в котором траектория встречается с xstar, как x-стоимость.или с вашим собственным методом интеграции.Я использовал odeint как документированный интерфейс, чтобы показать, как эта производная функция используется в интеграции.В результате время и значение у может быть экстремальным
ответил LutzL 14 J000000Saturday18 2018, 18:13:11

Похожие вопросы

Популярные теги

security × 330linux × 316macos × 2827 × 268performance × 244command-line × 241sql-server × 235joomla-3.x × 222java × 189c++ × 186windows × 180cisco × 168bash × 158c# × 142gmail × 139arduino-uno × 139javascript × 134ssh × 133seo × 132mysql × 132