模拟弹跳球?
Simulate a bouncing ball?
是否可以使用 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的时候进行干预,反转速度。还是我走错路了?
有点老套:
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).
来绘制弹跳球的轨迹
DifferentialEquations.jl offers sophisticated callbacks and event handling. Since the DifferentialEquations.jl algorithms are about 10x faster while offering a higher order interpolation,无论如何,这些算法显然是更好的选择。
第一个 link 是说明如何进行事件处理的文档。简单的界面使用宏。我从定义函数开始。
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
现在你照常解决。您使用 f
和初始条件定义 ODEProblem
,然后在时间跨度上调用 solve。唯一额外的是您将回调与求解器一起传递:
u0 = [50.0,0.0]
prob = ODEProblem(f,u0)
tspan = [0;15]
sol = solve(prob,tspan,callback=callback)
然后我们可以使用绘图配方自动绘制解决方案:
plot(sol)
结果是这样的:
这里有几点需要注意:
DifferentialEquations.jl 将自动使用插值来更安全地检查事件。例如,如果事件发生在一个时间步内但不在末尾,DifferentialEquations.jl 仍会找到它。可以将更多或更少的插值点作为选项包含在 @ode_event
宏中。
DifferentialEquations.jl 在事件发生的那一刻使用了寻根方法。即使自适应求解器跳过了事件,通过在插值上使用寻根,它也找到了事件的确切时间,从而正确处理了不连续性。您可以在图表中看到这一点,因为球永远不会变负。
它能做的还有很多。 Check out the docs。你几乎可以用它做任何事情。例如,让您的 ODE 改变 运行 模型中具有出生和死亡的细胞群的大小。这是其他求解器包做不到的。
即使具有所有这些功能,速度也不会受到影响。
如果您需要向 "ease of use" 界面宏添加任何额外功能,请告诉我。
是否可以使用 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的时候进行干预,反转速度。还是我走错路了?
有点老套:
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).
来绘制弹跳球的轨迹DifferentialEquations.jl offers sophisticated callbacks and event handling. Since the DifferentialEquations.jl algorithms are about 10x faster while offering a higher order interpolation,无论如何,这些算法显然是更好的选择。
第一个 link 是说明如何进行事件处理的文档。简单的界面使用宏。我从定义函数开始。
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
现在你照常解决。您使用 f
和初始条件定义 ODEProblem
,然后在时间跨度上调用 solve。唯一额外的是您将回调与求解器一起传递:
u0 = [50.0,0.0]
prob = ODEProblem(f,u0)
tspan = [0;15]
sol = solve(prob,tspan,callback=callback)
然后我们可以使用绘图配方自动绘制解决方案:
plot(sol)
结果是这样的:
这里有几点需要注意:
DifferentialEquations.jl 将自动使用插值来更安全地检查事件。例如,如果事件发生在一个时间步内但不在末尾,DifferentialEquations.jl 仍会找到它。可以将更多或更少的插值点作为选项包含在
@ode_event
宏中。DifferentialEquations.jl 在事件发生的那一刻使用了寻根方法。即使自适应求解器跳过了事件,通过在插值上使用寻根,它也找到了事件的确切时间,从而正确处理了不连续性。您可以在图表中看到这一点,因为球永远不会变负。
它能做的还有很多。 Check out the docs。你几乎可以用它做任何事情。例如,让您的 ODE 改变 运行 模型中具有出生和死亡的细胞群的大小。这是其他求解器包做不到的。
即使具有所有这些功能,速度也不会受到影响。
如果您需要向 "ease of use" 界面宏添加任何额外功能,请告诉我。