Вычисление повторного интеграла путем итеративного применения правила одномерных трапеций. ⇐ Python

[Anonymous]

У меня есть функция Python под названием catch_1D, которая вычисляет (приблизительный) интеграл одномерной функции f на интервале [a,b] с использованием правила составной трапеции. Теперь я хочу вычислить повторный интеграл в форме \int_{a}^{b} \int_{c}^{d} f(x,y) dy dx, итеративно вызываяrap_1D: Ниже представленrap_1D< /код>:

Код: Выделить всё

import numpy as np

def trap_1D(func, a, b, N, *args):

#Computes the integral of a univariate function over [a, b] using the composite trapezoidal rule.

h = (b - a) / N                #The grid spacing
x = np.linspace(a, b, N + 1)   #The sub-division points: x_0, x_1,...,x_N.
y = func(x, *args)

return (np.sum(y) - (y[0] + y[-1]) / 2) * h

Я протестировал catch_1D на нескольких функциях, и он работает. Далее я построил функцию I(x) := \int_{a}^{b} f(x,y) dy.

Код: Выделить всё

def create_I(f, a, b, N):
def I(x):
return trap_1D(lambda y: f(x, y), a, b, N)  #integrate the function y --> f(x,y)

return I

Я протестировал create_I, и он работает. Теперь, чтобы протестировать оценку повторных интегралов, я использую f(x,y) = sin(x) * sin(y) в качестве тестовой функции, которую я хочу проинтегрировать по прямоугольнику [0,pi/2] x [0 ,pi/2]:

Код: Выделить всё

def f(x,y):
return np.sin(x)*np.cos(y)

N = 100
I = create_I(f, 0.0, np.pi/2, N)   #The function I() is correct
I2 = trap_1D(I, 0.0, np.pi/2, N)   #This gives an error

Последняя строка выдает следующую ошибку:

Код: Выделить всё

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_15880\410820249.py in 
26 N = 100
27 I = create_I(f, 0.0, np.pi/2, 100)
---> 28 trap_1D(I, 0.0, np.pi/2, N)

~\AppData\Local\Temp\ipykernel_15880\410820249.py in trap_1D(func, a, b, N, *args)
9     y = func(x, *args)
10
---> 11     return (np.sum(y) - (y[0] + y[-1]) / 2) * h
12
13

IndexError: invalid index to scalar variable.

Я не понимаю эту ошибку. В чем здесь проблема?

Подробнее здесь: https://stackoverflow.com/questions/792 ... oidal-rule