如何在Python中求解积分方程？

Question

我正在尝试用Python解决这个积分方程：

其中z的范围是从0到1。

scipy.quad这个函数只能在某个区间内提供数值解，但它不能给出整个区间的解。

def f(z,Om,Ol): return 1./p.sqrt((1+z)**2 * (1+Om*z) - z*(2+z)*Ol)
quad(lambda r:f(r,Om,Ol),0,1)
(0.77142706642781111, 8.5645609096719596e-15)

不过我不知道怎么在这个区间内得到一个完整的向量，就像用scipy.odeint来解微分方程那样。

另一方面，sympy.integrate也做不到。我遇到了栈溢出的问题。而且，我也搞不清楚怎么把符号替换成一个列表，比如：

sy.integrate(x**2,x).subs(x,1)
1/3
sy.integrate(x**2,x).subs(x,[1,2])
TypeError: unhashable type: 'list'

所以我的问题是：有没有人知道怎么用Python解决一个积分方程？

Answer 1

在sympy的示例部分，http://docs.sympy.org/0.7.1/modules/integrals.html，他们展示了几乎相同问题的解决方案。你能把你的sympy代码发上来吗？

关于scipy，你有没有尝试用元组（tuple），它是可以被哈希的，而不是用列表（list）呢？比如：

sy.integrate(x**2,x).subs(x,(1,2,))

Answer 2

我想在这个积分前的 z 是个笔误，应该是 z1，你想要根据 DL 找到 z1。

首先，你得实现右边的部分 (rhs)：

def f(z,Om,Ol): 
        return 1./p.sqrt((1+z)**2 * (1+Om*z) - z*(2+z)*Ol)
def rhs(z1, Om, Ol, c, H0):
        return c/H0*(1+z1)*quad(lambda r:f(r, Om, Ol), 0, z1)[0]

接下来，你需要找到一个 z0，使得 rhs(z1, ...) = DL，这可以写成：

rhs(z1, ...) - DL = 0

这意味着你的问题可以简化为找到一个零点（只有一个，因为 rhs 是单调的），也就是：

f(z1) = rhs(z1, ...) - DL

在这里，你可以使用很多方法来找这个零点（可以参考 http://en.wikipedia.org/wiki/Root-finding_algorithm），还有 http://docs.scipy.org/doc/scipy/reference/optimize.html#root-finding 上的内容。

Answer 3

我明白你想要解决一个微分方程 dF/dz1 = f(z1, Om, Ol)，并且想在不同的位置得到 F(z1) 的值。如果是这样的话，SciPy 的 常微分方程（ODE）工具 是个不错的选择。你可能想要 看看 odeint()，因为它可以让你在指定的位置得到积分的值。