Регулирование вязкости с помощью решателей оды Сципи - PullRequest
1 голос
/ 21 апреля 2020

Для простоты рассмотрим следующее уравнение (уравнение Бюргерса):

Burgers equation

Давайте решим это с помощью scipy (в моем случае scipy.integrate.ode.set_integrator("zvode", ..).integrate(T)) с переменным решателем временных шагов.

Проблема заключается в следующем: , если мы используем наивную реализацию в пространстве Фурье

enter image description here

, тогда значение вязкости nu * d2x(u[t]) может вызвать превышение, если шаг по времени слишком велик. Это может привести к значительному шуму в решениях или даже к (фальшивым) расходящимся решениям (даже с жесткими решателями в несколько более сложной версии этого уравнения).

One способ для регуляризации этого состоит в том, чтобы оценить член вязкости на этапе t+dt, и этап обновления становится

enter image description here

Это решение хорошо работает, когда запрограммировано явно. Как я могу использовать переменный шаг Оде Солипа для его реализации? К моему удивлению, я не нашел никакой документации по этому довольно элементарному сложному вопросу ...

1 Ответ

2 голосов
/ 22 апреля 2020

Вы на самом деле не можете, или, с другой стороны, odeint или ode->zvode уже делает это с любой заданной проблемой.

Для первого вам нужно будет дать две части уравнение отдельно. Очевидно, что это не является частью интерфейса решателя. Посмотрите на решатели DDE и SDE, где такое разделение уравнения действительно требуется.

Во втором случае odeint и ode->zvode используют неявные многошаговые методы, что означает, что значения u(t+dt) и с правой стороны введите вычисление и лежащее в основе локальное приближение.

Вы все еще можете попытаться взломать свой исходный подход в решателе, предоставив функцию Якоби, которая содержит только второй производный член, но вполне вероятно, что вы это сделаете не достичь улучшения.

Вы можете разделить операторное ODE и решить линейную часть отдельно, введя

vhat(k,t) = exp(nu*k^2*t)*uhat(k,t)

так, чтобы

d/dt vhat(k,t) = -i*k*exp(nu*k^2*t)*conv(uhat(.,t),uhat(.,t))(k)
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...