2

我试图将这里给出的想法扩展为使用二维odeint的向量值函数的一维积分，其方法类似于.dblquad

下面你可以看到我目前的尝试：

import numpy as np
from scipy.integrate import odeint

def _infunc(x, func, gfun, hfun, more_args):
    a = gfun(x)
    b = hfun(x)
    y0 = f(x, a)
    return odeint(lambda v, y: f(x, y, *more_args), y0=y0, t=[a, b] )[1]

def dblodeint(f, a, b, gfun, hfun, args=()):
    y0 = f(a, gfun(a), *args)
    return odeint(lambda v, y: _infunc(y, f, gfun, hfun, args),
                  y0=y0, t=[a, b])[1]

if __name__ == '__main__':    
    def f(x, y):
        return np.array([x*y**2, x**2*y, x**4*y, x**6*y], float)

    def exact_int(a, b, ya, yb):
        return np.array([(b**2-a**2)*(yb**3-ya**3)/6.,
                         (b**3-a**3)*(yb**2-ya**2)/6.,
                         (b**5-a**5)*(yb**2-ya**2)/10.,
                         (b**7-a**7)*(yb**2-ya**2)/14.], float)

    print 'approx:', dblodeint(f, 0, 10, lambda x:0, lambda x:10)

    print 'exact:', exact_int(0, 10, 0, 10)

不幸的是，这不起作用......给出了以下错误的结果：

approx:Repeated convergence failures (perhaps bad Jacobian or tolerances).
Run with full_output = 1 to get quantitative information.
 [ 0.  0.  0.  0.]
exact: [  1.66666667e+04   1.66666667e+04   1.00000000e+06   7.14285714e+07]
4

0 回答 0