在朱莉娅中使用DifferentialEquations.jl处理不规则时间序列时，如何设置时间步长？

Question

利用谐振子的作用，在毫秒范围内由一个规则的时间序列w_i驱动微分方程。

ζ = 1/4pi                                           # damped ratio

function oscillator!(du,u,p,t)

      
           du[1] = u[2]                             # y'(t) = z(t)
           du[2] = -2*ζ*p(t)*u[2] - p(t)^2*u[1]     # z'(t) = -2ζw(t)z(t) -w(t)^2y(t)
       end
y0 = 0.0                                            # initial position
z0 = 0.0002                                         # initial speed
u0 = [y0, z0]                                       # initial state vector
tspan = (0.0,10)                                    # time interval

dt = 0.001                                          # timestep

w = t -> freq[Int(floor(t/dt))+1]                   # time series

prob = ODEProblem(oscillator!,u0,tspan,w)           # define ODEProblem

sol = solve(prob,DP5(),adaptive=false,dt=0.001)

当参数w_i是毫秒范围内的不规则时间序列时，如何设置时间步长。

date                    │  w  
────────────────────────┼───────
2022-09-26T00:00:00.023 │  4.3354
2022-09-26T00:00:00.125 │  2.34225
2022-09-26T00:00:00.383 │ -2.0312
2022-09-26T00:00:00.587 │ -0.280142 
2022-09-26T00:00:00.590 │  6.28319
2022-09-26T00:00:00.802 │  9.82271
2022-09-26T00:00:00.906 │ -5.21289
....................... |  ........

Answer 1

虽然可以禁用自适应，即使可以强制任意步长，但这通常不是您想要做的，因为它极大地限制了解决方案的准确性。

相反，将参数内插，使其接受任何t值。幸运的是，这样做真的很简单！

using Interpolations

...

ts = [0, 0.1, 0.4, 1.0]
ws = [1.0, 2.0, 3.0, 4.0]
w = linear_interpolation(ts, ws)
tspan = first(ts), last(ts)

prob = ODEProblem(oscillator!, u0, tspan, w)
sol = solve(prob, DP5(), dt=0.001)

当然，它不需要是线性插值。

如果仍然需要在特定时间点保存解决方案，请查看saveat for solve。例如，使用插值中使用的ts保存解决方案：

sol = solve(prob, DP5(), dt=0.001, saveat=ts)

编辑:跟进评论：

从数学上讲，你总是对整个领域的w(t)做一些假设。没有“由时间序列驱动”这样的说法。

例如，您在这里选择的标准Runge方法将要求ODE函数位于h/2。为了获得更好的DP5()，将在更多的子步骤中对其进行评估。这当然是不可避免的，不管是否使用适应性。尝试将println(t)添加到ODE函数中，您将看到以下内容。

如果有人来自matlab的ode45，而不是说它仍然简单地使用自适应，而只是像saveat那样处理显式的时间步骤。当然，它还将在显式步骤之外的各种t上评估函数。

因此，即使在第一个示例中，也是在插值w。您正在制造一种奇怪的constant_interpolation类型(但是使用floor，它与浮点数相结合会引起其他问题，因为floor(n*dt/dt)可能会评估为n或n-1)。

即使您要选择一种只在预定时间步骤(例如ExplicitEuler() )尝试计算的方法，您仍然在隐含地做出相同的假设，即w(t)在下一个时间步骤之前是常数。只是现在，您还从ODE集成中得到了一个更糟糕的解决方案。

如果前面的常量类型插值确实是在整个域中定义w(t)的方式(这就是您在floor(t/dt)中所做的)，那么我们所拥有的是：

w = extrapolate(interpolate((ts,), ws, Gridded(Constant{Previous}())), Flat())

在数学上，我们根本不可能忽略整个时间步骤发生的事情，也没有理由将时间限制在我们的"load“函数的样本点上。这在任何数学意义上都是不自然的。u'(t)必须在我们集成的整个域上定义。