解决弹道学的边界值问题
考虑一下奥运会上的铅球挑战。
- 给定:
-飞行时间
-初始位置 和
-结束位置 和
确定核心的初始速度和抛出的角度 和
弹道方程的推导
为了解决这个问题，有必要制定抛出物体的运动方程。:
根据该图，我们将绘制弹道运动方程。:
通过减少 ，并传递到一阶ODE系统，我们获得:
添加（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
const g = float(g_n).val
_,b,_,d = [1 2 3 4]
print("b=$b, d=$d")
让我们写下我们的方程组
作为函数 ballistic!(du,u,p,t).
p -系统参数, t -时间（但我们的系统是[自治]（https://ru.wikipedia.org/wiki/Автономная_система_дифференциальных_уравнений ）（静止），即它显然不依赖于 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
定义问题的边界条件
现在让我们介绍初始*(a)和最终(b)*时间点中的边界条件。
这些限制以以下形式介绍: .
也就是说，如果条件 然后
我们将根据铅球的世界纪录来获取任务值本身。:
- -基于第13页исследования
- -2023年奥运会纪录
- = 0
- （基于记录本身的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)
让我们对初始条件进行假设
α = 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]
求解两点边界值问题
让我们设置一个两点边界值问题（TwoPointBVProblem).
该问题由微分方程组、边界条件、初始假设和积分时间间隔组成。
关键词 bcresid_prototype 允许您设置初始假设的类型（我们可以在一个点上有4个条件 a 和4个条件在一个点 b. 但我们只指定了两个限制，因为速度在开始和结束时没有限制。
bvp = TwoPointBVProblem(ballistic!, (bca!, bcb!), u₀, tspan;
bcresid_prototype=(zeros(2), zeros(2)))
现在仍然是通过选择两种方法来解决问题 - Shooting，用求解器 Verns7(). 有关求解边界值问题的算法的详细信息，请参阅документации.
sol = solve(bvp, Shooting(Vern7()), saveat=0.1)
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系统的解确实满足边界条件。
现在我们找出投掷时刻速度矢量的模数和角度的值 .
norm([3 4]) # √(3²+4²)
norm([sol[1][3],sol[1][4]]) #√(v_x²+v_y²)
x |> f [相当于](https://engee.com/helpcenter/stable/ru/julia/base/base.html#Base.:|> ） f(x)
atan(sol[1][4]/sol[1][3]) |> rad2deg
因此，投掷的速度是≈ 14.5 м/с，以及投掷的角度？ 37.16°
确定投掷的角度，使您能够以最低的初始速度到达目标。
现在让我们来看看一组轨迹，并确保45°的角度允许您以最低的初始速度到达目标。
允许您根据初始速度更改轨迹颜色的功能
在构建轨迹之前，我们定义了函数，允许我们根据初始速度的值显示线的颜色。
我们来写一个函数 map_range，按比例显示段中的数字 [a,b] 段中 [A,B].
作为参数，它包含:
- 两个数字的元组
(a,b)-计划显示的初始部分 - 两个数字的元组
(A,B)-其中的部分
3计划显示。 -的数量[a,b]，计划在[A,B]
如果 ，然后函数抛出исключение。在我们的例子中，这个异常的类型是[DomainError](https://engee.com/helpcenter/stable/ru/julia/base/base.html#Core.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|
# map_range((14,18),(0,1),19) # раскомментируйте, чтобы проверить работу исключений
为了使用颜色和使用某些颜色集，必须连接库。 Colors 和 ColorSchemes
Pkg.add(["Colors","ColorSchemes"])
using Colors
using ColorSchemes
让我们创建一组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)
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），轨迹的自定义颜色变化