0

我求解了一组对 ODE，我使用与此示例类似的 GSL ODE 求解器求解。目前这是通过在python中编写文件来自动化的，例如

text = """
#include <stdio.h>
#include <gsl/gsl_errno.h>

...

"""

然后用相关单词替换字符串text并写入文件script.c。然后我用os.system它来运行它，例如

os.system("gcc -Wall -I/usr/local/include -c script.c")
os.system("gcc -static script.o -lgsl -lgslcblas -lm")
os.system("./a.out >  %s" % (file_name) )

所有这些都非常不雅，所以我一直在阅读有关替代品的文章，并且到目前为止偶然发现了 PyGSL、CythonGSL。这些似乎缺乏适当的文档，我还不够聪明，无法弄清楚它们是如何工作的！

任何建议将不胜感激。

4

2 回答 2

4

PyGSL 可能是一个合理的选择。（CythonGSL 也是，但我没试过。）

我正在运行 Ubuntu 12.04（64 位），对于这个安装，我使用的是 Anaconda python 发行版。（它可能适用于 Ubuntu 中的默认 Python，但我没有尝试。）

我从源代码安装了 GSL 1.15，并且我只是按照自述文件中的说明从源代码安装了 PyGSL 0.9.5。也就是说，我跑了

$ python setup.py build
$ python setup.py install

在提取的 tar 文件的顶层目录中。像许多构建 C 库的 python 包一样，在构建步骤期间会打印很多警告，但它完成时没有错误。

下面是一些演示在 PyGSL 中使用 ODE 求解器的代码：

#
# pendulum_demo.py
#
# Use PyGSL to solve the differential equations for a pendulum with
# friction.  Plot the solution using matplotlib.
#

import numpy as np
import matplotlib.pyplot as plt
from pygsl import odeiv


def pendsys(t, y, args):
    """
    The right-hand side of the first order system of differential equations.
    """
    theta, v = y
    g, b, L, m = args

    f = np.zeros((2,))
    f[0] = v
    f[1] = -b * v / (m * L ** 2) - np.sin(theta) * g / L

    return f


def pendjac(t, y, args):
    """
    The Jacobian of the system.
    """
    theta, v = y
    g, b, L, m = args

    jac = np.zeros((2, 2))
    jac[0, 1] = 1.0
    jac[1, 0] = -g / L * np.cos(theta)
    jac[1, 1] = -b / (m * L ** 2)

    dfdt = np.zeros((2,))

    return jac, dfdt


# Pendulum parameter values
#
# Gravitational constant
g = 9.81
# Friction coefficient
b = 0.5
# Pendulum length
L = 1.0
# Pendulum bob mass
m = 2.0

# Initial conditions
theta0 = np.pi - 0.01
v0 = 0.0
y = [theta0, v0]
t = 0

# Solver control parameters.
abserr = 1e-11
relerr = 1e-9

stoptime = 12.0

# Create the GSL ODE solver
N = 2
step    = odeiv.step_rk8pd(N, pendsys, pendjac, args=(g, b, L, m))
control = odeiv.control_y_new(step, abserr, relerr)
evolve  = odeiv.evolve(step, control, N)

# h is the initial step size.
h = stoptime / 500.0

# The time values and points in the solution are saved in the lists
# tsol and ysol, respectively.  The lists are initialized with
# the initial conditions.
tsol = [t]
ysol = [y]

# Call evolve.apply(...) until the solution reaches stoptime
while t < stoptime:
    t, h, y = evolve.apply(t, stoptime, h, y)
    tsol.append(t)
    ysol.append(y)

tsol = np.array(tsol)
ysol = np.array(ysol)

plt.plot(tsol, ysol[:, 0], label=r'$\theta$')
plt.plot(tsol, ysol[:, 1], label='v')
plt.xlabel('t')
plt.grid(True)
plt.legend()
plt.title('Pendulum with Friction')
plt.show()

结果图：

pygsl 生成的解图

于 2013-05-03T19:50:50.767 回答
1

If you can write a C library which takes parameters instead of using string substitution of C code from your Python code, then you can compile the C library seperately and call functions in it from your Python code using CFFI.

How complicated/dynamic are these substitutions?

于 2013-04-27T16:16:06.177 回答