python - 使用 python 内置函数进行耦合 ODE

Question

如果您需要的话，这部分只是背景

我正在开发二阶 Kuramoto 模型的数值求解器。下面给出了我用来求 theta 和 omega 导数的函数。

# 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。

# 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

感谢您的帮助!

Answer 1

您可以将现有函数包装到 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确实具有事件处理功能，但将其用于系统的事件收集相当麻烦。推进一些固定步骤，进行标准化和终止检测，然后重复此操作会更容易。

如果您想稍后提高效率，请直接使用 solve_ivp 后面的步进器类。