Hướng dẫn quad function in python - hàm quad trong python

Có gì sai khi chỉ tách nó ra thành phần thực và phần ảo? scipy.integrate.quadyêu cầu các phao trả về hàm tích hợp (còn gọi là số thực) cho thuật toán mà nó sử dụng.

import scipy
from scipy.integrate import quad

def complex_quadrature(func, a, b, **kwargs):
    def real_func(x):
        return scipy.real(func(x))
    def imag_func(x):
        return scipy.imag(func(x))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

Ví dụ,

>>> complex_quadrature(lambda x: (scipy.exp(1j*x)), 0,scipy.pi/2)
((0.99999999999999989+0.99999999999999989j),
 (1.1102230246251564e-14,),
 (1.1102230246251564e-14,))

đó là những gì bạn mong đợi để làm tròn lỗi - tích phân của exp (ix) từ 0, pi / 2 là (1 / i) (e ^ i pi / 2 - e ^ 0) = -i (i - 1) = 1 + i ~ (0,99999999999999989 + 0,99999999999999989j).

Và đối với hồ sơ trong trường hợp mọi người không rõ ràng 100%, tích phân là một hàm tuyến tính, có nghĩa là meaning {f (x) + kg (x)} dx = ∫ f (x) dx + k ∫ g (x ) dx (với k là hằng số đối với x). Hoặc đối với trường hợp cụ thể của chúng ta ∫ z (x) dx = ∫ Re z (x) dx + i ∫ Im z (x) dx as z (x) = Re z (x) + i Im z (x).

Nếu bạn đang cố gắng thực hiện tích hợp trên một đường dẫn trong mặt phẳng phức tạp (không phải dọc theo trục thực) hoặc vùng trong mặt phẳng phức tạp, bạn sẽ cần một thuật toán phức tạp hơn.

Lưu ý: Scipy.integrate sẽ không trực tiếp xử lý tích hợp phức tạp. Tại sao? Nó thực hiện công việc nặng nhọc trong thư viện FORTRAN QUADPACK , cụ thể là trong qagse.f , yêu cầu rõ ràng các hàm / biến phải có thực trước khi thực hiện "cầu phương thích ứng toàn cục dựa trên cầu phương Gauss – Kronrod 21 điểm trong mỗi subinterval, với sự gia tốc của Peter Thuật toán epsilon của Wynn. " Vì vậy, trừ khi bạn muốn thử và sửa đổi FORTRAN bên dưới để làm cho nó xử lý các số phức, biên dịch nó thành một thư viện mới, bạn sẽ không làm cho nó hoạt động.

Nếu bạn thực sự muốn thực hiện phương pháp Gauss-Kronrod với các số phức trong chính xác một phép tích phân, hãy xem trang wikipedias và thực hiện trực tiếp như thực hiện bên dưới (sử dụng quy tắc 15-pt, 7-pt). Lưu ý, hàm tôi memoize'd để lặp lại các lệnh gọi chung đến các biến chung (giả sử các lệnh gọi hàm chậm như thể hàm rất phức tạp). Cũng chỉ thực hiện quy tắc 7 pt và 15 pt, vì tôi không cảm thấy muốn tự mình tính toán các nút / trọng số và đó là những nút được liệt kê trên wikipedia, nhưng nhận được lỗi hợp lý cho các trường hợp thử nghiệm (~ 1e-14)

import scipy
from scipy import array

def quad_routine(func, a, b, x_list, w_list):
    c_1 = (b-a)/2.0
    c_2 = (b+a)/2.0
    eval_points = map(lambda x: c_1*x+c_2, x_list)
    func_evals = map(func, eval_points)
    return c_1 * sum(array(func_evals) * array(w_list))

def quad_gauss_7(func, a, b):
    x_gauss = [-0.949107912342759, -0.741531185599394, -0.405845151377397, 0, 0.405845151377397, 0.741531185599394, 0.949107912342759]
    w_gauss = array([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830050505119, 0.279705391489277,0.129484966168870])
    return quad_routine(func,a,b,x_gauss, w_gauss)

def quad_kronrod_15(func, a, b):
    x_kr = [-0.991455371120813,-0.949107912342759, -0.864864423359769, -0.741531185599394, -0.586087235467691,-0.405845151377397, -0.207784955007898, 0.0, 0.207784955007898,0.405845151377397, 0.586087235467691, 0.741531185599394, 0.864864423359769, 0.949107912342759, 0.991455371120813]
    w_kr = [0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525, 0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728, 0.204432940075298, 0.190350578064785, 0.169004726639267, 0.140653259715525,  0.104790010322250, 0.063092092629979, 0.022935322010529]
    return quad_routine(func,a,b,x_kr, w_kr)

class Memoize(object):
    def __init__(self, func):
        self.func = func
        self.eval_points = {}
    def __call__(self, *args):
        if args not in self.eval_points:
            self.eval_points[args] = self.func(*args)
        return self.eval_points[args]

def quad(func,a,b):
    ''' Output is the 15 point estimate; and the estimated error '''
    func = Memoize(func) #  Memoize function to skip repeated function calls.
    g7 = quad_gauss_7(func,a,b)
    k15 = quad_kronrod_15(func,a,b)
    # I don't have much faith in this error estimate taken from wikipedia
    # without incorporating how it should scale with changing limits
    return [k15, (200*scipy.absolute(g7-k15))**1.5]

Trường hợp thử nghiệm:

>>> quad(lambda x: scipy.exp(1j*x), 0,scipy.pi/2.0)
[(0.99999999999999711+0.99999999999999689j), 9.6120083407040365e-19]

Tôi không tin tưởng vào ước tính lỗi - Tôi đã lấy một cái gì đó từ wiki để ước tính lỗi được khuyến nghị khi tích hợp từ [-1 đến 1] và các giá trị có vẻ không hợp lý đối với tôi. Ví dụ: sai số ở trên so với sự thật là ~ 5e-15 chứ không phải ~ 1e-19. Tôi chắc rằng nếu ai đó tham khảo công thức nấu ăn num, bạn có thể nhận được ước tính chính xác hơn. (Có lẽ phải bội (a-b)/2số với một số quyền lực hoặc một cái gì đó tương tự).

Nhắc lại, phiên bản python kém chính xác hơn so với việc chỉ gọi tích hợp dựa trên QUADPACK của scipy hai lần. (Bạn có thể cải thiện nó nếu muốn).

54 hữu ích 4 bình luận chia sẻ 4 bình luận chia sẻ