如何从JiTCDDE实例获取导数值以调用add_past_points添加历史点

问题解答 #

不需要手动计算导数值，直接通过get_state()方法就能从JiTCDDE实例获取所有需要的状态和对应导数数据。

get_state()返回的序列每一项都是(时间戳, 状态数组, 导数数组)的三元组，和add_past_points()要求的输入参数完全匹配，你可以直接基于返回的数据做扰动修改后传入。

操作注意事项 #

• 必须在调用purge_past()之前执行get_state()获取所有历史锚点，清空后历史数据会丢失无法再获取。
• 你原有代码里的维度索引有误：前6维对应振荡器1应该取[:6]，你写的[6:]会拿到振荡器2的数值，注意不要搞反。
• 构造新锚点时，按你的扰动需求替换对应维度的状态和导数即可：要给振荡器1做10%周期的时间平移，就将每个时间点振荡器1的状态和导数，替换为该时间点减去偏移量时刻对应的值，振荡器2保持原有数值不变。
• 加完新锚点后需要调用adjust_diff()方法处理导数连续性，避免后续积分出现数值错误。
• 如果偏移量不是刚好对应已有锚点的时间戳，用线性插值处理状态和导数即可，稳态下插值误差可以忽略。

参考实现代码 #

import numpy as np

# 先完整获取当前所有历史锚点，必须在purge之前执行
old_anchors = DDE.get_state()
# 构建时间到锚点的映射，方便快速查询
time_list = [s[0] for s in old_anchors]
y_list = [s[1] for s in old_anchors]
dy_list = [s[2] for s in old_anchors]

new_anchors = []
for idx, (t, old_y, old_dy) in enumerate(old_anchors):
    # 计算要平移到的目标时刻
    target_t = t - perturbation
    # 插值获取目标时刻的状态和导数
    target_y = np.interp(target_t, time_list, y_list)
    target_dy = np.interp(target_t, time_list, dy_list)
    
    # 构造新的状态和导数：前6维用平移后的振荡器1数据，后6维保留原振荡器2数据
    new_y = np.hstack([target_y[:6], old_y[6:]])
    new_dy = np.hstack([target_dy[:6], old_dy[6:]])
    
    new_anchors.append( (t, new_y, new_dy) )

# 清空原有历史
DDE.purge_past()
# 加入扰动后的新锚点
DDE.add_past_points(new_anchors)
# 修正导数连续性
DDE.adjust_diff()
# 后续可正常调用integrate继续积分

内容的提问来源于stack exchange，提问作者Puff