对耦合的ODE使用python内置函数 [英] Using python built-in functions for coupled ODEs

问题描述

这部分只是您需要的背景

我正在为二阶仓本模型开发数值求解器.下面列出了我用于查找theta和omega导数的函数.

I am developing a numerical solver for the Second-Order Kuramoto Model. The functions I use to find the derivatives of theta and omega are given below.

# 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函数.

I know how to use odeint for two dimensions, (y,t). for my research I want to use a built-in Python function that works for higher dimensions.

问题2

我不一定要在预定时间后停止.我在下面的代码中还有其他停止条件，这些条件将指示方程组是否收敛到稳态.如何将它们整合到内置的Python解算器中?

I do not necessarily want to stop after a predetermined amount of time. I have other stopping conditions in the code below that will indicate whether the system of equations converges to the steady state. How do I incorporate these into a built-in Python solver?

我当前拥有的

这是我当前用于解决系统的代码.我刚刚在RK4中实现了恒定时间步进.

This is the code I am currently using to solve the system. I just implemented RK4 with constant time stepping in a loop.

# 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

感谢您的帮助！

推荐答案

您可以将现有函数包装到 odeint (选项 tfirst = True )接受的函数中，然后 solve_ivp 为

You can wrap your existing functions into a function accepted by odeint (option tfirst=True) and solve_ivp as

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函数.

The scipy methods pass numpy arrays and convert the return value into same, so that you do not have to care in the ODE function. The variant in comments is using explicit numpy functions.

虽然 solve_ivp 确实具有事件处理功能，但是将其用于系统的事件收集却相当麻烦.进行一些固定步骤，进行归一化和终止检测，然后重复此步骤会更容易.

While solve_ivp does have event handling, using it for a systematic collection of events is rather cumbersome. It would be easier to advance some fixed step, do the normalization and termination detection, and then repeat this.

如果以后要提高效率，请直接使用 solve_ivp 后面的步进器类.

If you want to later increase efficiency somewhat, use directly the stepper classes behind solve_ivp.

这篇关于对耦合的ODE使用python内置函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！