CODEWITHSUNDEEP
pythonscipysymbolic-mathnumerical-integrationnumerical-computing
Shankar_Dutt
Shankar_Dutt

Reputation: 125

How to evaluate double integral with variable upper limit of inner integral

I am trying to evaluate a double integral in which the limit of the inner integral is variable and depends on the outer integral. Coming from Matlab background, it is easy by using symbolic integration, but I have no idea how to calculate the same in Python in an efficient way. I read scipy documents but didn't get anything helpful. An example of such an integral is

Equation

Any help regarding the same would be very appreciated.

Update: I am looking for a numerical solution to the problem. Thanks a lot for the answers and help. They make complete sense. I have an updated question. I want to evaluate the following integral, which is an updated version of above. Thanks a lot in advance.

enter image description here

Upvotes: 1

Views: 1549

Answers (1)

Gorisanson
Gorisanson

Reputation: 2312

To compute a double integral, you can use the function scipy.integrate.dblquad. Or you can use the function scipy.integrate.quad twice.

The example integrals can be done like the following:

import math
from scipy.integrate import dblquad, quad
from scipy.special import erf, jv


def h(t, z):
    return f(t) * g(z)


def f(t):
    return 0.5 * t * (erf(t - a) - 1) * jv(0, q * t)


def g(z):
    return math.exp(-((z - a) ** 2)/(2 * (s ** 2)))


def h1(z):
    return integral_of_f(z) * g(z)


def h2(z):
    return (integral_of_f(z) ** 2) * g(z)


def integral_of_f(z):
    return quad(f, 0, 2 * z)[0]  # here abserr is discarded


if __name__ == '__main__':
    a, q, s = 0, 2, 3  # set the constants

    result, abserr = dblquad(h, 0, 60, lambda z: 0, lambda z: 2 * z)
    print(f'result: {result}, abserr: {abserr}')

    result1 = quad(h1, 0, 60)[0]  # abserr here for result1 is not valid since abserr is already discarded in the function integral_of_f
    print(f'result1: {result1}')

    result2 = quad(h2, 0, 60)[0]  # abserr here for result2 is not valid too for the same reason
    print(f'result2: {result2}')

which prints

result: -0.21705286423629183, abserr: 1.4361242750875712e-08
result1: -0.21705286423678177
result2: 0.013105371081178754

You can see that result, which is computed by using dblquad, and result1, which is computed by using quad twice, are almost the same value, because they computed the same integral.

result2 is for the updated question.

Upvotes: 4

Related Questions