使用 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

这些方程 return 个向量。

问题 1

我知道如何对二维 (y,t) 使用 odeint。对于我的研究，我想使用适用于更高维度的内置 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 数组并将return 值转换为相同的值，这样您就不必在ODE 函数中关心了。注释中的变体使用了显式的 numpy 函数。

虽然 solve_ivp 确实有事件处理功能，但使用它来系统地收集事件是相当麻烦的。提前一些固定的步骤，做归一化和终止检测，然后重复这个会更容易。

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

使用 python 耦合 ODE 的内置函数

Using python built-in functions for coupled ODEs

python

numpy

differential-equations

odeint