Цифровой двойной интеграция


Я сделал простую программу для численно aproximating двойной Интеграл, который соглашается с тем, что границы внутреннего интеграла-это функции:

import numpy as np
import time

def double_integral(func, limits, res=1000):
        t = time.clock()
        t1 = time.clock()
        t2 = time.clock()
        s = 0
        a, b = limits[0], limits[1]
        outer_values = np.linspace(a, b, res)
        c_is_func = callable(limits[2])
        d_is_func = callable(limits[3])
        for y in outer_values:
            if c_is_func:
                c = limits[2](y)
            else:
                c = limits[2]
            if d_is_func:
                d = limits[3](y)
            else:
                d = limits[3]
            dA = ((b - a) / res) * ((d - c) / res)
            inner_values = np.linspace(c, d, res)
            for x in inner_values:
                t2 = time.clock() - t2
                s += func(x, y) * dA
            t1 = time.clock() - t1
         t = time.clock() - t
         return s, t, t1 / res, t2 / res**2

Это, однако, ужасно медленно. Когда Рес=1000, таких, что интеграл-это сумма в миллион частей, это занимает около 5 секунд, чтобы добежать, но ответ будет только правильно, в 3-ем десятичное, по моему опыту. Есть ли способ ускорить этот процесс?

Код я использую, чтобы проверить Интеграл

def f(x, y):
    if (4 - y**2 - x**2) < 0:
        return 0              #This is to avoid taking the root of negarive #'s
    return np.sqrt(4 - y**2 - x**2)

def c(y):
    return np.sqrt(2 * y - y**2)

def d(y):
    return np.sqrt(4 - y**2)
# b d
# S S f(x,y) dx dy
# a c
a, b, = 0, 2
print(double_integral(f, [a, b, c, d]))

Интеграл eaqual в 16/9



Комментарии
3 ответа

Если вы хотите использовать numpyиспользуйте numpy правильно. Inestead из

for x in inner_values:
s += func(x, y) * dA

использовать более идиоматические, и гораздо быстрее

s += dA * np.sum(func(inner_values, y))

Примечание: это требует переписывания f как

return np.sqrt(np.maximum(0, 4 - y**2 - x**2))

поэтому он может принимать массив в качестве входных данных. Это не уменьшит точность, но приводит к гораздо более accaptable .04 секунд на 100-кратного улучшения. Основная идея заключается в пакете numpy-это не магия. Это обеспечивает быстрый векторизации.

6
ответ дан 22 февраля 2018 в 09:02 Источник Поделиться

Как отметил старик, в комментариях, двойной интеграции доступна в качестве составляющей scipy.integrate.dblquad. Это имеет аналогичный интерфейс код в пост:


scipy.integrate.dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

Вычислить двойной Интеграл.

Возврат двойной (определенный) интеграл func(y, x) от x = a..b и y = gfun(x)..hfun(x).

Различия (я) func берет свои аргументы в другом порядке; (II) в нижней и верхней границ должны быть указаны callables (но это не ограничение, потому что вы можете указать постоянные границы \$г=С\$ с функцией lambda x:c); (III) существует несколько аргументов для абсолютной и относительной толерантности результат.

Так для примера проблемы, которую вы хотели написать:

import numpy as np

def f(y, x):
return np.sqrt(np.maximum(0, 4 - y**2 - x**2))

def c(y):
return np.sqrt(2 * y - y**2)

def d(y):
return np.sqrt(4 - y**2)

а затем:

>>> import scipy.integrate
>>> scipy.integrate.dblquad(f, 0, 2, c, d)
(1.7777777777777706, 1.3374816809630374e-09)

(Второй результат является оценкой абсолютной погрешности в первом результате, несмотря на завышенную оценку в этом случае.)

1
ответ дан 23 февраля 2018 в 09:02 Источник Поделиться

Для функции strickly растет.
Между A и B
Рекурсивный поиск в точке Х interssection между линией A и B и сравнить Дельта Y, если слишком много вы добавить точку C и проверяем пересечение.
После этого вы должны суммировать трапеции sirfaces

-2
ответ дан 23 февраля 2018 в 11:02 Источник Поделиться