模拟一个弹跳球？

Question

是否有可能使用朱莉娅的方程求解器创建一个弹跳球的简单模型？

我从这开始:

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时进行干预,并反转速度.还是我走错了路？

Answer 1

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和初始条件,并在时间跨度上调用solve.唯一的补充是你将回调与求解器一起传递:

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

然后我们可以使用绘图配方自动绘制解决方案:

plot(sol)

结果是这样的:

这里有几点需要注意:

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

2. DifferentialEquations.jl使用rootfinding方法来磨练事件的时刻.即使自适应求解器逐步超过事件,通过在插值上使用rootfinding,它可以找到事件的确切时间,从而获得正确的不连续性.你可以在图表中看到,因为球永远不会消极.

3. 这可以做得更多.查看文档.你可以做任何事情.例如,让ODE在运行中改变大小,以模拟出生和死亡的细胞群.这是其他求解程序包无法做到的事情.

4. 即使具备所有这些功能,速度也不会受到影响.

如果您需要在"易用性"界面宏中添加任何额外功能,请告诉我们.