Каким будет вычислительно более быстрый способ реализовать это двумерное численное интегрирование?

Вы можете использовать Numba или низкоуровневую вызываемую

Почти ваш пример

Я просто передаю функцию непосредственно в scipy.integrate.dblquad вместо вашего метода, использующего лямбда-выражения для генерации функций.

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def f(t, z, q):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
        -0.5 * ((z - 40) / 2) ** 2)

def lower_inner(z):
    return 10.

def upper_inner(z):
    return 60.


y = np.empty(len(q))
for n in range(len(q)):
    y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()
print(end - start)
#143.73969149589539

Это уже немного быстрее (143 против 151 с), но единственная польза — простой пример для оптимизации.

Простая компиляция функций с помощью Numba

Чтобы это запустить, вам дополнительно понадобятся Numba и numba-scipy. Целью numba-scipy является предоставление обернутых функций из scipy.special.

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time
import numba as nb

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

#error_model="numpy" -> Don't check for division by zero
@nb.njit(error_model="numpy",fastmath=True)
def f(t, z, q):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
        -0.5 * ((z - 40) / 2) ** 2)

def lower_inner(z):
    return 10.

def upper_inner(z):
    return 60.


y = np.empty(len(q))
for n in range(len(q)):
    y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()
print(end - start)
#8.636585235595703

Использование вызываемого объекта низкого уровня

Функции scipy.integrate также предоставляют возможность передавать функцию обратного вызова C вместо функции Python. Эти функции можно написать, например, на C, Cython или Numba, которые я использую в этом примере. Основное преимущество заключается в том, что при вызове функции не требуется взаимодействие с интерпретатором Python.

Отличный ответ @Jacques Gaudin показывает простой способ сделать это, включая дополнительные аргументы.

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time
import numba as nb
from numba import cfunc
from numba.types import intc, CPointer, float64
from scipy import LowLevelCallable

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def jit_integrand_function(integrand_function):
    jitted_function = nb.njit(integrand_function, nopython=True)

    #error_model="numpy" -> Don't check for division by zero
    @cfunc(float64(intc, CPointer(float64)),error_model="numpy",fastmath=True)
    def wrapped(n, xx):
        ar = nb.carray(xx, n)
        return jitted_function(ar[0], ar[1], ar[2])
    return LowLevelCallable(wrapped.ctypes)

@jit_integrand_function
def f(t, z, q):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
        -0.5 * ((z - 40) / 2) ** 2)

def lower_inner(z):
    return 10.

def upper_inner(z):
    return 60.


y = np.empty(len(q))
for n in range(len(q)):
    y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()
print(end - start)
#3.2645838260650635

Как правило, выполнять суммирование с помощью матричных операций намного быстрее, чем использовать scipy.integrate.quad (или dblquad). Вы можете переписать свой f (q, z, t), чтобы взять вектор q, z и t и вернуть трехмерный массив f-значений, используя np.tensordot, а затем умножить свой элемент площади (dtdz) на значения функции и сумму их с помощью np.sum. Если ваш элемент области не является постоянным, вы должны создать массив элементов области и использовать np.einsum. Чтобы принять во внимание ваши пределы интеграции, вы можете использовать замаскированный массив для маскировки значений функции за пределами ваших пределов интеграции перед суммированием. Обратите внимание, что np.einsum игнорирует маски, поэтому, если вы используете einsum, вы можете использовать np.where, чтобы установить значения функций за пределами ваших ограничений интеграции равными нулю. Пример (с элементом постоянной площади и простыми пределами интегрирования):

import numpy as np
import scipy.special as ss
import time

def f(q, t, z):

    # Making 3D arrays before computation for readability. You can save some time by
    # Using tensordot directly when computing the output
    Mq = np.tensordot(q, np.ones((len(t), len(z))), axes=0)
    Mt = np.tensordot(np.ones(len(q)), np.tensordot(t, np.ones(len(z)), axes = 0), axes = 0)
    Mz = np.tensordot(np.ones((len(q), len(t))), z, axes = 0)

    return Mt * 0.5 * (ss.erf((Mt - Mz) / 3) - 1) * (Mq * Mt) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
     -0.5 * ((Mz - 40) / 2) ** 2)

q = np.linspace(0.03, 1, 1000)
t = np.linspace(0, 50, 250)
z = np.linspace(10, 60, 250)

#if you have constand dA you can shave some time by computing dA without using np.diff
#if dA is variable, you have to make an array of dA values and np.einsum instead of np.sum
t0 = time.process_time()
dA = np.diff(t)[0] * np.diff(z)[0]
func_vals = f(q, t, z)
I = np.sum(func_vals * dA, axis=(1, 2))
t1 = time.process_time()

это заняло 18,5 с на моем MacBook Pro 2012 года (2,5 ГГц i5) с dA = 0,04. Подобный подход также позволяет вам легко выбирать между точностью и эффективностью и устанавливать для dA значение, которое имеет смысл, когда вы знаете, как ведет себя ваша функция.

Тем не менее, стоит отметить, что если вы хотите получить большее количество очков, вы должны разделить свой интеграл, иначе вы рискуете максимально использовать свою память (1000 x 1000 x 1000) для удвоения требуется 8 ГБ оперативной памяти. Поэтому, если вы делаете очень большие интеграции с высокой точностью, может быть целесообразно быстро проверить требуемую память перед запуском.