Engee documentation
Notebook

Solution of the ballistics boundary value problem

Consider the problem of shot put in the Olympic Games.

  1. Danes:
  • flight time $(t_0,\:t_f)$
  • start positions $x(t_0)$ and $y(t_0)$
  • end positions $x(t_f)$ and $y(t_f)$

Determine the initial velocity of the cannonball and its throw angle $v_x(t_0)$ and $v_y(t_0)$

Derivation of the ballistics equation

To solve the problem, we need to derive the equations of motion of the thrown object: ballistics.png

Based on the figure, let's make the ballistics equations of motion: $$ \begin{split} &m \ddot{x} = 0\\ &m \ddot{y}=-mg \end{split} $$

Reducing $m$, and passing to the first order ODE system, we obtain: $$ \begin{split} &\dot{x} = v_x \\ &\dot{y} = v_y \\ &\dot{v}_x = 0 \\ &\dot{v}_y = -g \end{split} $$

Let's add (Pkg.add) and connect the libraries:

  • BoundaryValueDiffEq for solving the boundary value problem of differential equations
  • Plots - for displaying the graph of solutions
  • LinearAlgebra: norm to find the norm
  • PhysicalConstants will allow to use the value of standard acceleration of free fall
In [ ]: 
using Pkg
Pkg.add(["BoundaryValueDiffEq","LinearAlgebra","PhysicalConstants"])
using BoundaryValueDiffEq
using Plots
import LinearAlgebra: norm  # используем не весь модуль, а только определённую функцию
import PhysicalConstants.CODATA2018: g_n

g_n is not just a number, but a physical constant (value, dimensionality, uncertainty of the value). To get the "pure" value - you need to use float(g_n). It will return 9.80665 m s^-2. To get only the value, you need to access the val field.

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

Now let's look at technique, which allows you to insert "stubs"

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

b=2, d=4

Let's write our system of equations

$$ \begin{split} &\dot{x} = v_x \\ &\dot{y} = v_y \\ &\dot{v}_x = 0 \\ &\dot{v}_y = -g \end{split} $$ as a function of ballistic!(du,u,p,t). p - system parameters, t - time (but our system is autonomous (statsianar), i.e. it does not depend explicitly on t).

Don't be afraid to introduce new variables dx, ..., v_y for the sake of improving readability. Most of the time it doesn't affect performance.

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)

Defining the boundary conditions of the problem

Now let us introduce boundary conditions at the initial (a) and final (b) moments of time.

These boundary conditions are introduced in the following form: $g(u)=0$.

That is, if the condition $y(0)=2$, then $g(u)=y(0)-2=0$

The values of the problem will be taken based on the world record in shot put:

  • $x(0)=0$
  • $y(0)=2.1$ - based on page 13 research
  • $x(t_f)=23.1$ - Olympic record for 2023
  • $y(t_f)$ = 0
  • $(t_0,t_f)=(0,2)$ (based on the cropping of the shot_put.mp4 video of the record itself)
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)

Let's make assumptions about the initial conditions

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

Solution of the two-point boundary value problem

Let us define a two-point boundary value problem (TwoPointBVProblem).

This problem consists of a system of differential equations, boundary conditions, initial assumptions and integration time interval.

The keyword bcresid_prototype allows us to specify the kind of initial assumptions (we can have 4 conditions at a and 4 conditions at b. But we have specified only 2 constraints each, because the velocities at the beginning and the end were not constrained.

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

Now it remains to solve the problem, choosing as a method - Shooting, with the solver Verns7(). More information about algorithms of boundary value solvers can be found in documentation.

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]:

As can be seen from the figure, the solution of the ODE system does satisfy the boundary conditions.

Now let us find out the values of the modulus and angle of the velocity vector at the moment of throwing $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 equivalent f(x)

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

Thus, the throwing velocity ≈ 14.5 м/с, and the throwing angle ≈ 37.16°

Determine the throwing angle that allows you to reach the target with the lowest initial velocity.

Now let us consider a set of trajectories and verify that an angle of 45° allows us to reach the target with the lowest initial velocity.

Functions that allow you to change the colour of the trajectory, depending on the initial velocity

Before we plot the trajectories, let's define the functions that will allow us to display the colour of the line depending on the value of the initial velocity.

Let's write the function map_range, which proportionally maps the number from the segment [a,b] to the segment [A,B].

As parameters it contains:

  1. a tuple of two numbers (a,b) - the initial segment from which the mapping is to be performed
  2. tuple of two numbers (A,B) - the segment to which the mapping is planned.
  3. $x\in[a,b]$ - number from [a,b], which is planned to be mapped into the [A,B]

If $x\notin[a,b]$, the function throws exception. In our case, the type of this exception is DomainError, because x is out of the function definition area.

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) # раскомментируйте, чтобы проверить работу исключений

In order to work with colours and use certain sets of colours, it is necessary to connect libraries Colors and ColorSchemes

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

   Resolving package versions...
    Updating `~/.project/Project.toml`
⌅ [35d6a980] + ColorSchemes v3.26.0
⌅ [5ae59095] + Colors v0.12.11
  No Changes to `~/.project/Manifest.toml`

Let's create a set of 20 colours from blue to red:

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]]

Suppose we will have a set of trajectories whose initial velocities will vary from 14 to 17. And we want to find the colour corresponding to velocity 15, if blue is the lowest velocity and red is the highest.

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

Let's add empty graphs to which we will add new data in each iteration:

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_α)

Conclusions

We can conclude that indeed, in order to achieve the task at hand (throwing the ball a certain distance), an angle of 45° provides the minimum initial velocity at which this throw is possible. In other words, the 45° angle provides the maximum distance that the shot put can cover if the given initial velocity is fixed.

As a result of the task at hand, the following were considered:

  • solving the two-point boundary value problem
  • use of physical constants
  • creation of functions that throw exceptions when exceeding the boundaries of the function definition area
  • custom colour change of trajectories depending on the data corresponding to the trajectory (v_0, in our case)