我不明白如何使用 cython 设置和运行代码。
我在我的代码的相关部分添加了cdef,等。double
setup.py 当然没有使用名称 hello。 赛通文档
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
ext_modules = [Extension("hello", ["hello.pyx"])]
setup(
name = 'Hello world app',
cmdclass = {'build_ext': build_ext},
ext_modules = ext_modules
)
编辑
重新定义问题，使其更像一个具体的例子。
我想通过 Cython 运行这个 ODE 系统。我不确定是否double是最佳选择，但这些数字在 10 到 10 的数量级上，并且是非整数。
所以我想在更大的代码中调用它，其中mus在调用模块的脚本中定义了一个变量。我有我的setup.py文件来编译 pyx 文件。我不确定我需要做什么，所以我现在可以称之为颂歌。假设我将模块命名为 3bodyproblem。然后我会在 scipt 中调用它，import 3bodyproblem然后执行“3bodyproblem.3bodyproblem（这个输入是什么）”
我一直在阅读cython 的 odes 简介，但我不确定如何将他们的示例与我的示例一起使用。另外，如果需要是rk格式，请看第一个代码下面的代码。
代码 1
cdef deriv(double u, dt):
cdef double u[0], u[1], u[2]
return [u[3],
u[4],
u[5],
-mus * u[0] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
-mus * u[1] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
-mus * u[2] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5]
dt = np.linspace(0.0, t12sec, t12sec)
u = odeint(deriv, u0, dt, atol = 1e-13, rtol = 1e-13)
x, y, z, vx, vy, vz = u.T
代码 2
cdef deriv(dt, double u):
cdef double u[0], u[1], u[2], u[3], u[4], u[5]
return [u[3],
u[4],
u[5],
-mus * u[0] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
-mus * u[1] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
-mus * u[2] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5]
solver = ode(deriv).set_integrator('dop853')
solver.set_initial_value(u0)
z = np.zeros(300000)
for ii in range(300000):
z[ii] = solver.integrate(dt[ii])[0]
x, y, z, x2, y2, z2 = z.T
如果我只是尝试编译代码 1，cython 不知道定义较大.py文件的变量。我想避免在 cython 代码中设置变量，以便无需重新编译即可在其他地方使用它。以下是错误：
Error compiling Cython file:
------------------------------------------------------------
...
import numpy as np
from scipy.integrate import odeint
cdef deriv(double u, dt):
cdef double u[0], u[1], u[2]
^
------------------------------------------------------------
ODEcython.pyx:7:17: 'u' redeclared
Error compiling Cython file:
------------------------------------------------------------
...
import numpy as np
from scipy.integrate import odeint
cdef deriv(double u, dt):
cdef double u[0], u[1], u[2]
^
------------------------------------------------------------
ODEcython.pyx:7:23: 'u' redeclared
Error compiling Cython file:
------------------------------------------------------------
...
import numpy as np
from scipy.integrate import odeint
cdef deriv(double u, dt):
cdef double u[0], u[1], u[2]
^
------------------------------------------------------------
ODEcython.pyx:7:29: 'u' redeclared
Error compiling Cython file:
------------------------------------------------------------
...
-mus * u[0] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
-mus * u[1] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
-mus * u[2] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5]
dt = np.linspace(0.0, t12sec, t12sec)
^
------------------------------------------------------------
ODEcython.pyx:16:28: undeclared name not builtin: t12sec
Error compiling Cython file:
------------------------------------------------------------
...
-mus * u[1] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
-mus * u[2] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5]
dt = np.linspace(0.0, t12sec, t12sec)
u = odeint(deriv, u0, dt, atol = 1e-13, rtol = 1e-13)
^
------------------------------------------------------------
ODEcython.pyx:17:16: Cannot convert 'object (double, object)' to Python object
Error compiling Cython file:
------------------------------------------------------------
...
-mus * u[1] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
-mus * u[2] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5]
dt = np.linspace(0.0, t12sec, t12sec)
u = odeint(deriv, u0, dt, atol = 1e-13, rtol = 1e-13)
^
------------------------------------------------------------
ODEcython.pyx:17:20: undeclared name not builtin: u0
Error compiling Cython file:
------------------------------------------------------------
...
cdef deriv(double u, dt):
cdef double u[0], u[1], u[2]
return [u[3],
u[4],
u[5],
-mus * u[0] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
^
------------------------------------------------------------
ODEcython.pyx:11:17: undeclared name not builtin: mus
building 'ODEcython' extension
creating build
creating build/temp.linux-x86_64-2.7
x86_64-linux-gnu-gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC -I/usr/include/python2.7 -c ODEcython.c -o build/temp.linux-x86_64-2.7/ODEcython.o
ODEcython.c:1:2: error: #error Do not use this file, it is the result of a failed Cython compilation.
error: command 'x86_64-linux-gnu-gcc' failed with exit status 1
也许这会有所帮助
所以我想在这样的文件中使用模块：（不同的派生函数，但忽略这个想法是一样的）
import numpy as np
from scipy.integrate import ode
import pylab
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import brentq
from scipy.optimize import fsolve
me = 5.974 * 10 ** 24 # mass of the earth
mm = 7.348 * 10 ** 22 # mass of the moon
G = 6.67259 * 10 ** -20 # gravitational parameter
re = 6378.0 # radius of the earth in km
rm = 1737.0 # radius of the moon in km
r12 = 384400.0 # distance between the CoM of the earth and moon
rs = 66100.0 # distance to the moon SOI
Lambda = np.pi / 6 # angle at arrival to SOI
M = me + mm
d = 300 # distance the spacecraft is above the Earth
pi1 = me / M
pi2 = mm / M
mue = 398600.0 # gravitational parameter of earth km^3/sec^2
mum = G * mm # grav param of the moon
mu = mue + mum
omega = np.sqrt(mu / r12 ** 3)
# distance from the earth to Lambda on the SOI
r1 = np.sqrt(r12 ** 2 + rs ** 2 - 2 * r12 * rs * np.cos(Lambda))
vbo = 10.85 # velocity at burnout
h = (re + d) * vbo # angular momentum
energy = vbo ** 2 / 2 - mue / (re + d) # energy
v1 = np.sqrt(2.0 * (energy + mue / r1)) # refer to the close up of moon diagram
# refer to diagram for angles
theta1 = np.arccos(h / (r1 * v1))
phi1 = np.arcsin(rs * np.sin(Lambda) / r1)
#
p = h ** 2 / mue # semi-latus rectum
a = -mue / (2 * energy) # semi-major axis
eccen = np.sqrt(1 - p / a) # eccentricity
nu0 = 0
nu1 = np.arccos((p - r1) / (eccen * r1))
# Solving for the eccentric anomaly
def f(E0):
return np.tan(E0 / 2) - np.sqrt((1 - eccen) / (1 + eccen)) * np.tan(0)
E0 = brentq(f, 0, 5)
def g(E1):
return np.tan(E1 / 2) - np.sqrt((1 - eccen) / (1 + eccen)) * np.tan(nu1 / 2)
E1 = fsolve(g, 0)
# Time of flight from r0 to SOI
deltat = (np.sqrt(a ** 3 / mue) * (E1 - eccen * np.sin(E1)
- (E0 - eccen * np.sin(E0))))
# Solve for the initial phase angle
def s(phi0):
return phi0 + deltat * 2 * np.pi / (27.32 * 86400) + phi1 - nu1
phi0 = fsolve(s, 0)
nu = -phi0
gamma = 0 * np.pi / 180 # angle in radians of the flight path
vx = vbo * (np.sin(gamma) * np.cos(nu) - np.cos(gamma) * np.sin(nu))
# velocity of the bo in the x direction
vy = vbo * (np.sin(gamma) * np.sin(nu) + np.cos(gamma) * np.cos(nu))
# velocity of the bo in the y direction
xrel = (re + 300.0) * np.cos(nu) - pi2 * r12
yrel = (re + 300.0) * np.sin(nu)
u0 = [xrel, yrel, 0, vx, vy, 0]
def deriv(u, dt):
return [u[3], # dotu[0] = u[3]
u[4], # dotu[1] = u[4]
u[5], # dotu[2] = u[5]
(2 * omega * u[4] + omega ** 2 * u[0] - mue * (u[0] + pi2 * r12) /
np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum *
(u[0] - pi1 * r12) /
np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
# dotu[3] = that
(-2 * omega * u[3] + omega ** 2 * u[1] - mue * u[1] /
np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum * u[1] /
np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
# dotu[4] = that
0] # dotu[5] = 0
dt = np.linspace(0.0, 259200.0, 259200.0) # secs to run the simulation
u = odeint(deriv, u0, dt)
x, y, z, x2, y2, z2 = u.T
fig = pylab.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z, color = 'r')
# adding the moon
phi = np.linspace(0, 2 * np.pi, 100)
theta = np.linspace(0, np.pi, 100)
xm = rm * np.outer(np.cos(phi), np.sin(theta)) + r12 - pi2 * r12
ym = rm * np.outer(np.sin(phi), np.sin(theta))
zm = rm * np.outer(np.ones(np.size(phi)), np.cos(theta))
ax.plot_surface(xm, ym, zm, color = '#696969', linewidth = 0)
ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])
# adding the earth
xe = re * np.outer(np.cos(phi), np.sin(theta)) - pi2 * r12
ye = re * np.outer(np.sin(phi), np.sin(theta))
ze = re * np.outer(np.ones(np.size(phi)), np.cos(theta))
ax.plot_surface(xe, ye, ze, color = '#4169E1', linewidth = 0)
ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])
pylab.show()