0

我正在使用 JiTCDDE 模块来求解延迟微分方程。我的工作要求我让模型发展一段时间，然后对其进行扰动。为此，我尝试使用jitcdd.purge_past()后跟jitcdde.add_past_points(). 问题是，后者要求我不仅提供输入时的值，还需要提供导数的值。

有没有办法从jitcdde实例本身获取这些，还是我需要手动笨拙地计算它们？

编辑：更多信息

我的系统由两个耦合在一起的非线性振荡器组成（它们不是相位振荡器）。我让系统发展一段时间，直到它达到稳定状态，然后通过稍微移动两个振荡器中的一个来扰乱它。I shift it 的量计算为振荡周期的百分比，即有效相移。

我的系统是 12 维的，我正在努力获得一个最小的工作示例，所以这是代码的非工作模型：

f = [...]
DDE = jitcdde(f)
DDE.constant_past([1]*len(f))
DDE.step_on_discontinuities()

initial_integration_time = 30
DDE.integrate(initial_integration_time)

这样做之后，我想执行应该是的扰动，比如说一个周期的 10%，假设我知道静止周期长度是T。我要做的是用来get_state获取系统和导数的当前状态，以及状态和导数perturbation时间单位之前的状态。然后，我可以构造一组新的anchors，我可以将其传递给它DDE.add_past_points以从那里开始集成。

T = ...
perturbation = 0.1*T #in units of time

state = DDE.get_state()
# get position and derivative of first oscilator
# since the system is 12-D, the first six correspond to osc1
osc1_x = [state[1][6:] for s in state]
osc1_dx = [state[2][6:] for s in state]

# here, do the same for osc2 at DDE.t - perturbation
4

1 回答 1

1

问题

要回答发布的问题，要从 jitcdde 实例中获取导数，您需要做的就是调用该get_state方法。假设DDE您的jitcdde实例已经集成：

state = DDE.get_state()
derivatives = [s[2] for s in state]

这是有效的，因为get_state将返回一个Cubic Hermite Spline实例，该实例的行为基本上类似于一个元组列表（加上一些很酷的方法），其中每个元组包含一个时间、当时系统的状态以及当时系统的导数，请参阅文档中的此条目以获取更多信息。

我的问题

为了具体解决我的问题，如果任何路人碰巧关心，我做了以下（与问题相同的模型样式的代码）：

# define and initialize system
f = [...]
DDE = jitcdde(f)
DDE.constant_past([1]*len(f))
DDE.step_on_discontinuities()

initial_integration_time = 30
T = ... # period length
perturbation = 0.1*T # 10% of period, in units of time

# get state of system at perturbation time units before the 
# end of the initial integration
DDE.integrate(initial_integration_time-perturbation)
state1 = DDE.get_state.copy()

# get state after initial integration is done
DDE.integrate(initial_integration_time)
state2 = DDE.get_state()

# generate new anchors from these set of states
eps = 1e-6
anchors = []
for i, s in enumerate(state2[::-1]): #iterate in reverse
    # stop if I reached the end, there's no telling 
    # which states is gonna be longer
    if i >= len(state1):
        break
    
    # calculate new anchor at the time of the anchor I'm looking at
    # shifted by the perturbation length
    tt = s[0] - perturbation_duration
    state1.interpolate_anchor(tt)
    _, x1, dx1 = state1[state1.last_index_before(tt+eps)]
    
    x_new = np.hstack( (x1[:6], s[1][6:]) )
    dx_new = np.hstack( (dx1[:6], s[2][6:]) )
    anchors.append( (tt, x_new, dx_new) )

# add this new past to the integrator and evolve the system
DDE.purge_past()
DDE.add_past_points(anchors)
DDE.step_on_discontinuities()

# now I can evolve the system as I please

这段代码有一些注意事项，比如确保扰动不会太大并且状态足够长，但这是我的解决方案背后的基本思想。我不会详细解释它的作用，因为我认为它对任何人都没有指导意义。

于 2021-11-04T21:29:08.793 回答