2

是否可以使用 Julia 的方程求解器创建一个简单的弹跳球模型？

我从这个开始：

using ODE

function bb(t, f)
    (y, v) = f
    dy_dt = v
    dv_dt = -9.81
    [dy_dt, dv_dt]
end

const y0 =  50.0             # height
const v0 =   0.0             # velocity
const startpos = [y0; v0]

ts = 0.0:0.25:10             # time span

t, res = ode45(bb, startpos, ts)

这会产生有用的数字：

julia> t
44-element Array{Float64,1}:
  0.0
  0.0551392
  0.25
  0.5
  0.75
  1.0
  ⋮
  8.75
  9.0
  9.25
  9.5
  9.75
 10.0

julia> res
44-element Array{Array{Float64,1},1}:
 [50.0,0.0]
 [49.9851,-0.540915]
 [49.6934,-2.4525]
 [48.7738,-4.905]
 [47.2409,-7.3575]
 ⋮
 [-392.676,-93.195]
 [-416.282,-95.6475]
 [-440.5,-98.1]

但不知何故，它需要在高度为 0 时进行干预，并反转速度。还是我走错了路？

4

2 回答 2

5

DifferentialEquations.jl 提供复杂的回调和事件处理。由于DifferentialEquations.jl 算法在提供更高阶插值的同时快了大约 10 倍，因此这些算法显然是更好的选择。

第一个链接是显示如何进行事件处理的文档。简单的界面使用宏。我从定义函数开始。

f = @ode_def BallBounce begin
  dy =  v
  dv = -g
end g=9.81

在这里，我展示了ParameterizedFunctions.jl以使语法更好，但您可以将函数直接定义为就地更新f(t,u,du)（如 Sundials.jl）。接下来，您定义确定事件何时发生的函数。它可以是任何正数并在事件时间达到零的函数。在这里，我们正在检查球何时落地，或者何时y=0，所以：

function event_f(t,u) # Event when event_f(t,u,k) == 0
  u[1]
end

接下来你说当事件发生时要做什么。在这里，我们想要反转速度的符号：

function apply_event!(u,cache)
  u[2] = -u[2]
end

您将这些函数放在一起以使用宏构建回调：

callback = @ode_callback begin
  @ode_event event_f apply_event!
end

现在你照常解决。您定义ODEProblem使用f条件和初始条件，并在时间跨度上调用求解。唯一额外的是您将回调与求解器一起传递：

u0 = [50.0,0.0]
prob = ODEProblem(f,u0)
tspan = [0;15]
sol = solve(prob,tspan,callback=callback)

然后我们可以使用 plot recipe 自动绘制解决方案：

plot(sol)

结果是这样的：

弹跳

这里有几点需要注意：

  1. DifferentialEquations.jl 将自动使用插值来更安全地检查事件。例如，如果事件发生在一个时间步内但不是在结束时，DifferentialEquations.jl 仍会找到它。可以包含更多或更少的插值点作为@ode_event宏的选项。

  2. DifferentialEquations.jl 使用寻根方法来磨练事件发生的时刻。即使自适应求解器越过事件，通过在插值上使用求根，它也可以找到事件的确切时间，从而获得正确的不连续性。您可以在图表中看到这一点，因为球永远不会变成负数。

  3. 这可以做的还有很多。查看文档。你可以用这个做几乎任何事情。例如，让您的 ODE 在运行过程中改变大小，以模拟具有出生和死亡的细胞群。这是其他求解器包无法做到的。

  4. 即使拥有所有这些功能，速度也不会受到影响。

如果您需要在“易用性”界面宏中添加任何额外功能，请告诉我。

于 2016-10-16T01:56:45.910 回答
3

有点hacky：

function bb(t, f)
    (y, v) = f
    dy_dt = v
    dv_dt = -9.81*sign(y)
    [dy_dt, dv_dt]
end

您只需遵循约定，其中 y 和 -y 指的是相同的高度。然后，您可以通过绘制 abs(y) 来绘制弹跳球的轨迹。

于 2016-10-14T01:10:29.900 回答