弹道边界值问题的求解¶

考虑奥林匹克运动会中的铅球问题。

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)$

弹道方程的推导¶

为了解决这个问题，我们需要推导出被抛物体的运动方程：

根据图示，我们来建立弹道运动方程： $$ \begin{split} &m \ddot{x} = 0\\ &m \ddot{y}=-mg \end{split} $$

还原$m$ ，并转入一阶 ODE 系统，得到 $$ \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 允许使用标准自由落体加速度值

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 字段。

g_n

Standard acceleration of gravitation (g_n)
Value                         = 9.80665 m s^-2
Standard uncertainty          = (exact)
Relative standard uncertainty = (exact)
Reference                     = CODATA 2018

const g = float(g_n).val

9.80665

现在让我们来看看技术，它允许您插入 "存根"。

_,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 - 时间（但我们的系统是自主的 (statsianar)，即它并不明确依赖于t ）。

不要害怕为了提高可读性而引入新变量dx, ..., v_y 。大多数情况下，这不会影响性能。

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

ballistic! (generic function with 1 method)

定义问题的边界条件¶

现在，让我们引入初始 (a) 时刻和最终 (b) 时刻的边界条件。

这些边界条件的引入形式如下：$g(u)=0$.

也就是说，如果条件为$y(0)=2$ ，那么$g(u)=y(0)-2=0$

问题的取值将以铅球世界纪录为基础：

-$x(0)=0$ -$y(0)=2.1$ - 基于第 13 页[研究]（https://www.researchgate.net/publication/326479265_Biomechanical_Report_for_the_IAAF_World_Championships_2017_Shot_Put_Men's） -$x(t_f)=23.1$ - 2023 年奥运会纪录 -$y(t_f)$ = 0 -$(t_0,t_f)=(0,2)$ （基于对记录本身的 shot_put.mp4 视频的裁剪）

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)

(0.0, 2)

让我们对初始条件进行假设

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

4-element Vector{Float64}:
 0.0
 2.1
 7.0710678118654755
 7.071067811865475

两点边界值问题的解法¶

让我们定义一个两点边界值问题 (TwoPointBVProblem) 。

这个问题由微分方程系、边界条件、初始假设和积分时间间隔组成。

关键字bcresid_prototype 允许我们指定初始假设的种类（在a 可以有 4 个条件，在b 可以有 4 个条件。但我们只分别指定了 2 个约束条件，因为开始和结束时的速度不受约束。

bvp = TwoPointBVProblem(ballistic!, (bca!, bcb!), u₀, tspan;
        bcresid_prototype=(zeros(2), zeros(2)))

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

Verns7()现在只剩下用Shooting 的求解器求解这个问题了。有关边界值求解器算法的更多信息，请参阅文档。

sol = solve(bvp, Shooting(Vern7()), saveat=0.1)

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]

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="краевые условия")

从图中可以看出，ODE 系统的解确实满足边界条件。

现在让我们求出抛出时刻速度矢量的模量和角度值$t=0$ 。

norm([3 4]) # √(3²+4²)

5.0

norm([sol[1][3],sol[1][4]])  #√(v_x²+v_y²)

14.49418570401595

x |> f 等效f(x)

atan(sol[1][4]/sol[1][3]) |> rad2deg

37.16764031481598

因此，投掷速度 ≈14.5 м/с ，投掷角度 ≈37.16°

确定能以最小初速度到达目标的投掷角度。¶

现在，让我们考虑一组轨迹，并验证 45° 角能使我们以最小的初速度到达目标。

根据初速度改变轨迹颜色的函数¶

在绘制轨迹之前，我们先来定义一些函数，这些函数可以根据初速度值来显示轨迹线的颜色。

我们先来写一个函数map_range ，它可以按比例将数字从线段[a,b] 映射到线段[A,B] 。

参数包括

两个数字的二进制数(a,b) - 要进行映射的初始线段
两个数字的元组(A,B) - 计划进行映射的分段。
$x\in[a,b]$ - 来自[a,b] 的数字，计划将其映射到[A,B]

如果$x\notin[a,b]$ ，函数将抛出异常。在本例中，该异常的类型是 DomainError，因为x 不在函数定义区域内。

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|

0.25

# map_range((14,18),(0,1),19) # раскомментируйте, чтобы проверить работу исключений

为了处理颜色和使用某些颜色集，有必要连接Colors 和ColorSchemes

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`

让我们创建一组从蓝色到红色的 20 种颜色：

COLORS_NUMBER = 20
ndr = ColorScheme(get(ColorSchemes.darkrainbow, range(0.0, 1.0, length=COLORS_NUMBER)))

[ndr[1],ndr[end]]

假设我们有一组轨迹，其初始速度从 14 到 17。如果蓝色是最低速度，红色是最高速度，那么我们要找到与速度 15 相对应的颜色。

round(Int,3.5)

4

ndr[round(Int,map_range((14,17),(1,COLORS_NUMBER),15))]

让我们添加空图表，每次迭代时都添加新数据：

plt_trajectories = plot()
scatter_v_0_α = scatter();

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)

display(scatter_v_0_α)

结论¶

我们可以得出这样的结论：为了完成手头的任务（将球投掷到一定距离），45° 角提供了实现投掷的最小初速度。换句话说，如果给定的初速度固定不变，45° 角提供了推铅球所能达到的最大距离。

根据手头的任务，我们考虑了以下几点：

解决两点边界值问题
使用物理常数
创建在超出函数定义区域边界时抛出异常的函数
根据轨迹对应的数据（本例中为 v_0）自定义轨迹的颜色变化