使用 Cython 的 RK4 代码运行速度并不比平常的 python 快

问题描述 投票:0回答:0

我想编写一个简单的代码来使用龙格-库塔方法求解耦合微分方程。用python写的代码太长了,所以我用Cython写的,就是

import time
import cython
cimport numpy as np

def f(np.ndarray[np.float64_t, ndim=1] state, float t):
    cdef float x, p
    cdef np.ndarray[np.float64_t, ndim=1] du
    x = state[0]; p = state[1];
    du = np.array([p,-x])
    return du

def motion(np.ndarray[np.float64_t, ndim=1] u0, np.ndarray[np.float64_t, ndim=1] t):
    cdef float dt, tmax, t0
    dt = t[1]-t[0]; tmax = t[-1]; t0 = t[0]
    cdef int nt = int((tmax-t0)/(dt));
    cdef np.ndarray[np.float64_t, ndim=1] Xt = np.array((nt+1)*[list(u0)[0]]), Pt = np.array((nt+1)*[list(u0)[1]]), u
    cdef np.ndarray[np.float64_t, ndim=1] k1,k2,k3,k4

    u = u0
    for ii in range(0,nt):
        #print(i)
        k1 = dt * f(u, t[ii])
        k2 = dt * f(u + 0.5 * k1, t[ii] + 0.5*dt)
        k3 = dt * f(u + 0.5 * k2, t[ii] + 0.5*dt)
        k4 = dt * f(u + k3, t[ii] + dt)
        #u[ii+1] = u[ii] + (k1 + 2*(k2 + k3) + k4)/6.0

        Xt[ii+1] = Xt[ii] + (k1 + 2*(k2 + k3) + k4)[0]/6.0
        Pt[ii+1] = Pt[ii] + (k1 + 2*(k2 + k3) + k4)[1]/6.0

        u = np.array( [ Xt[ii+1],Pt[ii+1] ] )

    return Xt,Pt

此代码可以正确编译,但是编译所需的时间与原始 python 代码相同。如何更正代码?

python cython
© www.soinside.com 2019 - 2024. All rights reserved.