Использование python встроенных функций для связанных ODE - PullRequest
Новая
       80

Использование python встроенных функций для связанных ODE

0 голосов
/ 24 февраля 2020

ЭТА ЧАСТЬ ПРОСТО ФОН, ЕСЛИ ВАМ НУЖНО ЭТО

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

# n-dimensional change in omega
def d_theta(omega):
    return omega

# n-dimensional change in omega
def d_omega(K,A,P,alpha,mask,n):
    def layer1(theta,omega):
        T = theta[:,None] - theta
        A[mask] = K[mask] * np.sin(T[mask])
        return - alpha*omega + P - A.sum(1)
    return layer1

Эти уравнения возвращают векторы.

ВОПРОС 1

Я знаю, как использовать odeint для двух измерений (y, t). для своих исследований я хочу использовать встроенную функцию Python, которая работает для более высоких измерений.

ВОПРОС 2

Не обязательно останавливаться по истечении заданного промежутка времени. У меня есть другие условия остановки в коде ниже, который укажет, сходится ли система уравнений к установившемуся состоянию. Как мне включить их во встроенный Python решатель?

ЧТО Я В ТЕЧЕНИЕ В настоящее время

Это код, который я в настоящее время использую для решения системы. Я только что внедрил RK4 с постоянным шагом по времени в al oop. 

# This function randomly samples initial values in the domain and returns whether the solution converged

# Inputs:
# f             change in theta (d_theta)
# g             change in omega (d_omega)
# tol           when step size is lower than tolerance, the solution is said to converge
# h             size of the time step
# max_iter      maximum number of steps Runge-Kutta will perform before giving up
# max_laps      maximum number of laps the solution can do before giving up
# fixed_t       vector of fixed points of theta
# fixed_o       vector of fixed points of omega
# n             number of dimensions

# theta         initial theta vector
# omega         initial omega vector

# Outputs:
# converges     true if it nodes restabilizes, false otherwise

def kuramoto_rk4_wss(f,g,tol_ss,tol_step,h,max_iter,max_laps,fixed_o,fixed_t,n):
    def layer1(theta,omega):
        lap = np.zeros(n, dtype = int)
        converges = False
        i = 0
        tau = 2 * np.pi

        while(i < max_iter): # perform RK4 with constant time step
            p_omega = omega
            p_theta = theta
            T1 = h*f(omega)
            O1 = h*g(theta,omega)
            T2 = h*f(omega + O1/2)
            O2 = h*g(theta + T1/2,omega + O1/2)
            T3 = h*f(omega + O2/2)
            O3 = h*g(theta + T2/2,omega + O2/2)
            T4 = h*f(omega + O3)
            O4 = h*g(theta + T3,omega + O3)

            theta = theta + (T1 + 2*T2 + 2*T3 + T4)/6 # take theta time step

            mask2 = np.array(np.where(np.logical_or(theta > tau, theta < 0))) # find which nodes left [0, 2pi]
            lap[mask2] = lap[mask2] + 1 # increment the mask
            theta[mask2] = np.mod(theta[mask2], tau) # take the modulus

            omega = omega + (O1 + 2*O2 + 2*O3 + O4)/6

            if(max_laps in lap): # if any generator rotates this many times it probably won't converge
                break
            elif(np.any(omega > 12)): # if any of the generators is rotating this fast, it probably won't converge
                break
            elif(np.linalg.norm(omega) < tol_ss and # assert the nodes are sufficiently close to the equilibrium
               np.linalg.norm(omega - p_omega) < tol_step and # assert change in omega is small
               np.linalg.norm(theta - p_theta) < tol_step): # assert change in theta is small
                converges = True
                break
            i = i + 1
        return converges
    return layer1

Спасибо за вашу помощь!

1 Ответ

1 голос
/ 24 февраля 2020

Вы можете превратить ваши существующие функции в функцию, принятую odeint (опция tfirst=True) и solve_ivp как 

def odesys(t,u):
    theta,omega = u[:n],u[n:]; # or = u.reshape(2,-1);
    return [ *f(omega), *g(theta,omega) ]; # or np.concatenate([f(omega), g(theta,omega)])

u0 = [*theta0, *omega0]
t = linspan(t0, tf, timesteps+1);
u = odeint(odesys, u0, t, tfirst=True);

#or

res = solve_ivp(odesys, [t0,tf], u0, t_eval=t)

Методы scipy передают numpy массивы и преобразуют верните значение в то же самое, так что вам не нужно заботиться о функции ODE. Вариант в комментариях использует явные функции numpy.

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

Если позже вы захотите немного повысить эффективность, используйте непосредственно классы степпера за solve_ivp.

...