Реализация Matlab «OutputFcn» в Python для решения ODE - PullRequest
Новая
       16

Реализация Matlab «OutputFcn» в Python для решения ODE

0 голосов
/ 28 марта 2019

Я пытаюсь решить ODE, используя Python solve_ivp.Однако я хочу динамически изменить правую часть моего ODE на основе сравнения между текущим решением и предыдущим решением.Идея этого заключается в том, что моя правая часть - это векторное поле, и я хочу обеспечить направленность векторного поля, перевернув правую часть на основе направления предыдущего решения.

Реализация для этого заключается в следующем: я хочу проверить скалярное произведение в определении функции справа между предыдущим решением и векторным полем.Если скалярное произведение отрицательно, правая часть умножается на -1.

Поэтому мне нужно получить доступ к предыдущему состоянию решателя ODE и использовать его по сравнению с текущей итерацией.В MATLAB есть возможность использования «OutputFcn» при решении ODE.Эта функция вызывается после каждой итерации интегратора.Поэтому в функции можно просто извлечь состояние как переменную и использовать его на следующей итерации.Мне не удалось найти что-то подобное для Python. 

def RHS(timesnotused,x):
    out = solve_ivp(doubleGyreVar, [0,T/2, T], [x[0], x[1], 1, 0, 0, 1], rtol = 1e-10, atol=1e-10)
    output = out.y
    J = output[2:,-1].reshape(2,2)
    CG = np.matmul(J.T , J)
    lambdas, xis = np.linalg.eig(CG)
    xi_1 = xis[np.argmin(lambdas)] 
    xi_2 = xis[np.argmax(lambdas)]
    lambda_1 = np.min(lambdas)
    lambda_2 = np.max(lambdas)

    alpha = ((lambda_2-lambda_1) / (lambda_2+lambda_1))**2

    sign = 1
    if np.dot(xlast,xi_1) < 0:
        sign = -1

    return(sign*alpha*xi_1)

Как видно, я хочу, чтобы "xlast" было предыдущим решением, и проверяю его с помощью xi_1 текущей итерации.Каким-то образом xlast необходимо обновлять каждую итерацию.

1 Ответ

0 голосов
/ 28 марта 2019

Если дифференциальное уравнение, которое вы хотите решить, является 

dx/dt = sign(g(x)) * F(x)

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

Самый ясный метод для решения такой многофазной системы состоит в том, чтобы представить только непрерывные правые части числовому решателю и обработать изменение фазы с помощью механизма событий, который также присутствует в scipy.integrate.solve_ivp.

Я исследовал механизм, чтобы сделать это с помощью инструментов scipy в аналогичной проблеме со знаковой функцией, создающей разрыв в odeint возвращает неправильные результаты для ODE, включая функцию дискретного

...