Почему scipy.odeint с разным размером шага в массиве t дает качественно разные результаты?
Я исследую странные аттракторы и в настоящее время работаю над аттрактором Даффинга (странный 2D-аттрактор). Траектория решения делает петли в положительном и отрицательном направлении x, и меня особенно интересует последовательность, в которой она переключается между ними (паттерн хаотичен). Это довольно легко оценить, поскольку вы можете построить график зависимости значения x траектории от времени (независимая переменная). Я использую следующий код:
a=0.35
b=0.3
w=1
def dX_dt(X, t):
return np.array([X[1],
X[0]-(X[0]**3)-a*X[1]+b*np.cos(w*t)])
tstepmax=5000
tmax=1000
tmin=0
t1 = np.linspace(tmin, tmax+tmin, tstepmax)
X1=odeint(dX_dt,[0,0],t1)
Я не был уверен, сколько шагов требуется для заданного количества единиц времени, поэтому экспериментировал с настройкойtstepmax
до 50 000 и 500000. К моему удивлению, через довольно короткое время они дали качественно разные результаты, как видно на этом графике:Три наложенные траектории. У меня сложилось впечатление, чтоt1
вектор был просто местом, где траектория оценивалась для целей построения, и метод внутренне определял размер шага, который нужно использовать; действительно, для некоторых методов размер шага не является постоянным на траектории. Однако это изменение дает принципиально иную траекторию. Почему? И, что, возможно, более важно, какова «настоящая» траектория? Как я могу достоверно получить «настоящую» траекторию?
1 ответ
Аналогичная проблема наблюдалась и обсуждалась в https://scicomp.stackexchange.com/questions/42613/need-help-to-full-understand-scipys-odeints-reported-step-sizes-eval-times.
Короче говоря, причина, по-видимому, в том, что вызов степпера odeint (в коде lsoda на Фортране) не выполняет полный шаг, а только до следующего выходного значения. Состояние степпера сохраняется в массиве рабочей области, но, видимо, не полностью. Некоторые управляющие переменные, кажется, сбрасываются при каждом вызове, каждый вызов фактически создает новый шаговый объект, смешивая старую информацию и информацию по умолчанию. Это приводит к несколько разной внутренней последовательности шагов и, следовательно, к несколько разной эволюции ошибок, что в хаотической системе может привести к заметно расходящимся решениям.
С использованиемsolve_ivp
вам не следует наблюдать эту проблему, поскольку временной цикл выполняет полные внутренние шаги, а затем интерполирует все соответствующие выходные значения для текущего интервала. Класс Stepper продолжает существовать со всей внутренней информацией. Таким образом, внутренняя последовательность шагов должна быть независимой отt_eval
аргумент.