Nicholaswogan / numbalsoda

Python wrapper of LSODA (solving ODEs) which can be called from within numba functions.
MIT License
91 stars 9 forks source link

Can I use 2d-array of 'u0'? #18

Closed kimyoungjin06 closed 1 year ago

kimyoungjin06 commented 1 year ago

I want to use the 2d-shape initial values of the multiple-particles system. But, I think I cannot use the slicing or indexing.

def swing_lsoda(network):
    """
    rhs = swing_cover(network)
    funcptr = rhs.address
    ...
    data = np.array([1.0])
    usol, success = lsoda(funcptr, u0, t_eval, data = data)
    """
    @nb.cfunc(lsoda_sig) # network
    def rhs(t, u, du, p): # p0=m, p1=gamma, p2=P, p3=K
        # u is 2d-array 
        Interaction = p[3] * SinIntE(network, u[0])
        du[0] = u[1] #Thetas of nodes
        du[1] = 1/p[0]*(p[2] - p[1]*u[1] - Interaction) #Omega 

How can I try it?

Nicholaswogan commented 1 year ago

u0 can not be a 2d array. But you can flatten the 2D array, then reshape it in the rhs function. For example

@cfunc(lsoda_sig)
def rhs(t, u, du, p):
    u_2D = nb.carray(u, (1,2))
    # rest of function goes here
    du[0] = u_2D[0,0]-u_2D[0,0]*u_2D[0,1]
    du[1] = u_2D[0,0]*u_2D[0,1]-u_2D[0,1]*p[0]

funcptr = rhs.address
u0_2d = np.array([[5.,0.8]])
u0 = u0_2d.flatten()
data = np.array([1.0])
t_eval = np.linspace(0.0,50.0,1000)

usol, success = lsoda(funcptr, u0, t_eval, data = data)
kimyoungjin06 commented 1 year ago

@Nicholaswogan Perfect! ;)

Thanks a lot! I think this is so beginner level in the numba, but I also put in the first step in the numba.

kimyoungjin06 commented 1 year ago

@Nicholaswogan

I have an additional question.

Can I not use slicing?

def swing_lsoda(network, N):
    rhs = swing_cover(network)
    funcptr = rhs.address
    ...
    data = np.array([1.0])
    usol, success = lsoda(funcptr, u0, t_eval, data = data)

    @nb.cfunc(lsoda_sig) # network
    def rhs(t, u, du, p): # p0=m, p1=gamma, p2=P, p3=K
        # u_2D = nb.carray(u, (N,2))
        Interaction = p[3] * SinIntE(network, u[:N])
        du[0] = u[N:] #Thetas of nodes
        du[1] = 1/p[0]*(p[2] - p[1]*u[N:] - Interaction) #Omegas of nodes
Nicholaswogan commented 1 year ago

You can slice, but only after you have converted u to a numpy array with nb.carray

@nb.cfunc(lsoda_sig) # network
def rhs(t, u, du, p): 
    u_2D = nb.carray(u,(N,2))
    # `u_2D` can be sliced, but `u` can not be
    u_1D = nb.carray(u,(N*2))
    # `u_1D` can be sliced as well