我想编写一个简单的代码来使用龙格-库塔方法求解耦合微分方程。用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 代码相同。如何更正代码?