Документация Engee
Notebook

Решение краевой задачи баллистики

Открыть пример в Engee

Рассмотрим задачу толкания ядра на олимпийских играх.

  1. Даны:
  • время полёта $(t_0,\:t_f)$
  • начальные положения $x(t_0)$ и $y(t_0)$
  • конечные положения $x(t_f)$ и $y(t_f)$

Определить начальную скорость ядра и угол его броска $v_x(t_0)$ и $v_y(t_0)$

Вывод уравнения баллистики

Для решения задачи необходимо составить уравнения движения брошенного предмета: ballistics.png

Исходя из рисунка, составим уравнения движения баллистики: $$ \begin{split} &m \ddot{x} = 0\\ &m \ddot{y}=-mg \end{split} $$

Сократив $m$, и перейдя к системе ОДУ первого порядка, получим: $$ \begin{split} &\dot{x} = v_x \\ &\dot{y} = v_y \\ &\dot{v}_x = 0 \\ &\dot{v}_y = -g \end{split} $$

Добавим (Pkg.add) и подключим библиотеки:

  • BoundaryValueDiffEq для решения краевой задачи дифференциальных уравнений
  • Plots - для отображения графика решений
  • LinearAlgebra: norm для нахождения нормы
  • PhysicalConstants поззволит использовать значение стандартного ускорения свободного падения
In [ ]: 
using Pkg
Pkg.add(["BoundaryValueDiffEq","LinearAlgebra","PhysicalConstants"])
using BoundaryValueDiffEq
using Plots
import LinearAlgebra: norm  # используем не весь модуль, а только определённую функцию
import PhysicalConstants.CODATA2018: g_n

g_n представляет собой не просто число, а физическую константу (значение, размерность, неопределённость величины). Чтобы получить "чистое" значение - нужно воспользоваться float(g_n). Она вернёт 9.80665 m s^-2. Чтобы получить только значение, нужно обратиться к полю val

In [ ]: 
g_n
Out[0]:
Standard acceleration of gravitation (g_n)
Value                         = 9.80665 m s^-2
Standard uncertainty          = (exact)
Relative standard uncertainty = (exact)
Reference                     = CODATA 2018
In [ ]: 
const g = float(g_n).val
Out[0]:
9.80665

Теперь рассмотрим приём, который позволяет вставлять "заглушки"

In [ ]: 
_,b,_,d = [1 2 3 4]
print("b=$b, d=$d")

b=2, d=4

Запишем нашу систему уравнений

$$ \begin{split} &\dot{x} = v_x \\ &\dot{y} = v_y \\ &\dot{v}_x = 0 \\ &\dot{v}_y = -g \end{split} $$ в виде функции ballistic!(du,u,p,t). p - параметры системы, t - время (но наша система является автономной (стацианарной), т.е. явно не зависит от t).

Не бойтесь вводить новые переменные dx, ..., v_y ради повышения читабельности. Чаще всего это не влияет на производительность.

In [ ]: 
function ballistic!(du, u, p, t)
    _, _, v_x, v_y = u
    
    du[1] =  dx  = v_x
    du[2] =  dy  = v_y
    du[3] = dv_x = 0
    du[4] = dv_y = -g
end
Out[0]:
ballistic! (generic function with 1 method)

Определение краевых условий задачи

Теперь введём граничные условия (boundary conditions) в начальный (a) и конечный (b) моменты времени.

Эти ограничения вводятся в следующем виде: $g(u)=0$.

Т.е. если условие $y(0)=2$, то $g(u)=y(0)-2=0$

Сами значения задачи будем брать, исходя из мирового рекорда по толканию ядра:

  • $x(0)=0$
  • $y(0)=2.1$ - исходя из стр 13 исследования
  • $x(t_f)=23.1$ - олимпийский рекорд на 2023 год
  • $y(t_f)$ = 0
  • $(t_0,t_f)=(0,2)$ (исходя из кадрирования видео shot_put.mp4 самого рекорда)
In [ ]: 
function bca!(resid_a, u_a, p)
    resid_a[1] = u_a[1]        #   x(0) = 0
    resid_a[2] = u_a[2] - 2.1  #   y(0) = 2
end

function bcb!(resid_b, u_b, p)
    resid_b[1] = u_b[1] - 23.1  # x(t_f) = 95
    resid_b[2] = u_b[2]       # y(t_f) = 0
end

tspan = (0., 2)
Out[0]:
(0.0, 2)

Сделаем предположения о начальных условиях

In [ ]: 
α = deg2rad(45)
x₀ = 0
y₀ = 2.1
v₀ = 10
v_x_0 = v₀ * cos(α)
v_y_0 = v₀ * sin(α)
u₀ = [x₀, y₀, v_x_0, v_y_0]
Out[0]:
4-element Vector{Float64}:
 0.0
 2.1
 7.0710678118654755
 7.071067811865475

Решение двухточечной краевой задачи

Зададим двухточечную краевую задачу (TwoPointBVProblem).

Эта проблема состоит из системы дифференциальных уравнений, краевых условий, начальных предположений и временного интервала интегрирования.

Ключевое слово bcresid_prototype позволяет задать вид начальных предположений (у нас может быть 4 условия в точке a и 4 условия в точке b. Но мы указали только по 2 ограничения, так как скорости и в начале и конце не ограничивали.

In [ ]: 
bvp = TwoPointBVProblem(ballistic!, (bca!, bcb!), u₀, tspan;
        bcresid_prototype=(zeros(2), zeros(2)))
Out[0]:
BVProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 2.0)
u0: 4-element Vector{Float64}:
 0.0
 2.1
 7.0710678118654755
 7.071067811865475

Теперь осталось решить задачу, выбрав в качетсве метода - Shooting, с решателем Verns7(). Подробнее про алгоритмы решателей краевой задачи можно прочитать в документации.

In [ ]: 
sol = solve(bvp, Shooting(Vern7()), saveat=0.1)
Out[0]:
retcode: Success
Interpolation: 1st order linear
t: 21-element Vector{Float64}:
 0.0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1.0
 1.1
 1.2
 1.3
 1.4
 1.5
 1.6
 1.7
 1.8
 1.9
 2.0
u: 21-element Vector{Vector{Float64}}:
 [0.0, 2.1, 11.55, 8.756650000000016]
 [1.1549999999999925, 2.926631749999996, 11.55, 7.7759850000000235]
 [2.3100000000000005, 3.6551970000000034, 11.55, 6.795320000000016]
 [3.465, 4.285695750000003, 11.55, 5.814655000000018]
 [4.619999999999997, 4.818128000000005, 11.55, 4.8339900000000195]
 [5.774999999999997, 5.252493750000007, 11.55, 3.85332500000002]
 [6.9300000000000015, 5.588793000000004, 11.55, 2.872660000000015]
 [8.084999999999988, 5.827025749999997, 11.55, 1.8919950000000263]
 [9.24000000000005, 5.967191999999978, 11.55, 0.9113299999999747]
 [10.395000000000012, 6.009291750000039, 11.55, -0.0693349999999926]
 [11.549999999999898, 5.953325000000021, 11.55, -1.049999999999896]
 [12.704999999999917, 5.799291749999795, 11.55, -2.030664999999912]
 [13.860000000000117, 5.547191999999763, 11.55, -3.011330000000079]
 [15.014999999999933, 5.197025750000018, 11.55, -3.9919949999999247]
 [16.170000000000005, 4.7487929999997505, 11.55, -4.972659999999985]
 [17.325000000000557, 4.202493749999437, 11.55, -5.953325000000454]
 [18.4800000000015, 3.5581279999992574, 11.55, -6.933990000001257]
 [19.635000000000005, 2.815695750000004, 11.55, -7.914654999999986]
 [20.790000000000003, 1.9751970000000036, 11.55, -8.895319999999986]
 [21.945, 1.0366317500000068, 11.55, -9.875984999999982]
 [23.100000000000005, 5.879140493979664e-15, 11.55, -10.856649999999986]
In [ ]: 
gr()  # используем статическое отображение графиков
plt_record = plot(  sol[1, :], sol[2, :],
                    aspect_ratio=1.0, xlabel="x, м",ylabel="y, м", label = "траектория")
scatter!(plt_record, [0,23.1],[2.1, 0],label="краевые условия")
Out[0]:

Как видно из рисунка, решение системы ОДУ действительно удовлетворяет краевым условиям.

Теперь узнаем значения модуля и угла вектора скорости в момент броска $t=0$.

In [ ]: 
norm([3 4]) # √(3²+4²)
Out[0]:
5.0
In [ ]: 
norm([sol[1][3],sol[1][4]])  #√(v_x²+v_y²)
Out[0]:
14.49418570401595

x |> f эквивалентно f(x)

In [ ]: 
atan(sol[1][4]/sol[1][3]) |> rad2deg
Out[0]:
37.16764031481598

Таким образом, скорость броска ≈ 14.5 м/с, и угол броска ≈ 37.16°

Определение угла броска, который позволяет достигать цели при наименьшей начальной скорости.

Теперь рассмотрим набор траекторий, и убедимся в том, что угол в 45° позволяет достигнуть цели при наименьшей начальной скорости.

Функции, позволяющие изменять цвет траектории, в зависимости от начальной скорости

Перед тем как построить траектории, определим функции, которые позволят отображать цвет линии в зависимости от значения начальной скорости.

Напишем функцию map_range, пропорционально отображающую число из отрезка [a,b] в отрезок [A,B].

В качестве параметров она содержит:

  1. кортеж из двух чисел (a,b) - исходный отрезок, из которого планируется отображение
  2. кортеж из двух чисел (A,B) - отрезок, в который планируется отображение
  3. $x\in[a,b]$ - число из [a,b], которое планируется отобразить в [A,B]

Если $x\notin[a,b]$, то функция выбрасывает исключение. В нашем случае тип этого исключения - DomainError, так как x выходит из области определения функции.

In [ ]: 
map_range((a, b), (A, B), x) = a < x < b ?
                           A + (x - a) / (b - a) * (B - A) :
                           throw(DomainError(x, "x must be in domain ($a, $b)"))
map_range((14,18),(0,1),15)

# |14---15___16___17___18| |0---0.25___0.5___0.75___1|
Out[0]:
0.25
In [ ]: 
# map_range((14,18),(0,1),19) # раскомментируйте, чтобы проверить работу исключений

Для того чтобы работать с цветами и использовать определённые наборы цветов необходимо подключить библиотеки Colors и ColorSchemes

In [ ]: 
Pkg.add(["Colors","ColorSchemes"])
using Colors
using ColorSchemes

Создадим набор из 20 цветов, расположеных от синего к красному:

In [ ]: 
COLORS_NUMBER = 20
ndr = ColorScheme(get(ColorSchemes.darkrainbow, range(0.0, 1.0, length=COLORS_NUMBER)))
Out[0]:
No description has been provided for this image
In [ ]: 
[ndr[1],ndr[end]]
Out[0]:
No description has been provided for this image

Предположим, что у нас будет набор траекторий, начальные скорости которых будут меняться от 14 до 17. И мы хотим найти цвет, соответствующий скорости 15, если синий - наименьшая скорость, красный - наибольшая.

In [ ]: 
round(Int,3.5)
Out[0]:
4
In [ ]: 
ndr[round(Int,map_range((14,17),(1,COLORS_NUMBER),15))]
Out[0]:
No description has been provided for this image

Добавим пустые графики, на которые будем в каждой итерации добавлять новые данные:

In [ ]: 
plt_trajectories = plot()
scatter_v_0_α = scatter();
In [ ]: 
for tspan in map(x -> (0, x), 1.5:0.15:3)
    # решение краевой задачи для заданного (t_0,t_f)
    bvp = TwoPointBVProblem(ballistic!, (bca!, bcb!), u₀, tspan;
        bcresid_prototype=(zeros(2), zeros(2)))
    sol = solve(bvp, Shooting(Vern7()), saveat=0.1)
    # определение v_0 и α 
    v_0 = round(norm([sol[1][3],sol[1][4]]),digits=2)   # (round, π, digits=4) = 3.1416
    α = atan(sol[1][4]/sol[1][3]) |> rad2deg
    # построение графиков
    ## построение траекторий
    plot!(  plt_trajectories, sol[1, :], sol[2, :], aspect_ratio=1.0,
            color=ndr[round(Int,map_range((14,17), (1,COLORS_NUMBER),v_0))],
            label=nothing, w=3,xlabel="x",ylabel="y", title="траектории толкания ядра")  
    ## для построения зависимости v_0(α)
    scatter!(scatter_v_0_α,[α],[v_0],label=nothing,
             xlabel="α", ylabel="v_0",title="Зависимость начальной скорости от угла толкания",
              color=ndr[round(Int,map_range((14,17), (1,COLORS_NUMBER),v_0))])  
end
display(plt_trajectories)
In [ ]: 
display(scatter_v_0_α)

Выводы

Можем сделать вывод, что действительно, для достижения поставленной задачи (метания ядра на определённое расстояние) угол в 45° обеспечивает минимальную начальную скорость, при которой этот бросок возможен. Иными словами, угол в 45° обеспечивает максимальное расстояние, которое может преодолеть ядро, если заданная начальная скорость фиксирована.

В результате выполения поставленной задачи были рассмотрены:

  • решение двухточечной краевой задачи
  • использование физических констант
  • создание функций, выбрасывающих исключения при выходе за границы области определения функции
  • пользовательское изменение цвета траекторий в зависимости от данных, соответствующих траектории (v_0, в нашем случае)