如何在python中求解积分方程？答案

【问题标题】：How to solve an integral equation in python?如何在python中求解积分方程？
【发布时间】：2011-11-01 22:05:30
【问题描述】：

我正在尝试使用 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 求解积分方程吗？

【问题讨论】：

使用 numpy：docs.scipy.org/doc/numpy/reference/generated/…
当您谈论积分方程时：您要寻找的值是多少？你想得到给定 D_L 的 z1 吗？
积分前有一个 z，这是错字吗？
我不知道如何使用 linalg 求解积分方程。

标签： python integration scipy equation

【解决方案1】：

我了解到您想要求解一个微分方程 dF/dz1 = f(z1, Om, Ol) 并希望 F(z1) 在不同的位置。如果是这种情况，那么 SciPy 的Ordinary Differential Equation (ODE) routines 就是要走的路。您可能特别想checkodeint()，因为它可以在您指定的位置为您提供积分值。

【讨论】：

感谢您的帮助，但我不是要解微分方程，而是要解积分方程。我在其他情况下使用了 odeint() 例程，它对我来说很好，但它只对微分方程有用，在这种情况下我不能使用它。
@Illa Rivero Losada：谢谢。您是否想让 Python 找到积分的分析公式？ “求解积分方程”实际上具有非常具体的含义（en.wikipedia.org/wiki/Integral_equation），这似乎不适用于您的问题。如果您确实在寻找分析解决方案，那么提出这个问题的好地方就是 mathoverflow。 ΩM = 0 的情况是可行的（en.wikipedia.org/wiki/List_of_integrals_of_irrational_functions），但我不确定一般情况......

【解决方案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://docs.scipy.org/doc/scipy/reference/optimize.html#root-finding 的多种寻根方法（参见http://en.wikipedia.org/wiki/Root-finding_algorithm）

【讨论】：

【解决方案3】：

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

对于 scipy，您是否尝试过使用可散列的元组而不是列表？例如：

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

【讨论】：

我试过了，但出现错误：TypeError: unsupported operand type(s) for ** or pow(): 'tuple' and 'Integer'

【解决方案4】：

我终于使用了带有 for 语句的“quad”函数来解决我的问题：

import pylab as p
import scipy as s
from scipy.integrate import odeint,quad

def Dl_lflrw(z,Om,Ol):
    c = 1.
    H0 = 1.
    y = []
    for i in z:
        y1 = (c/H0)*(1+i)*quad(lambda r:f(r,Om,Ol),0,i)[0]
        y.append(y1)
    return p.asarray(y)

def f(z,Om,Ol):
    return 1./p.sqrt((1+z)**2 * (1+Om*z) - z*(2+z)*Ol)

感谢大家的想法，他们真的很有帮助。

【讨论】：