通过迭代应用一维梯形规则来计算迭代积分

Question

我有一个名为

trap_1D

的 Python 函数，它使用复合梯形规则计算单变量函数 f 在区间 [a,b] 上的（近似）积分。现在我想通过迭代调用

trap_1D

来计算 \int_{a}^{b} \int_{c}^{d} f(x,y) dy dx 形式的迭代积分：下面是

trap_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

我已经在一些功能上测试了

trap_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 <module>
     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.

我不明白这个错误。是什么导致了这里的问题？

Answer 1

您的

create_I

函数应该返回 I 的向量化版本。目前，当 linspace x 传递到函数中时，当您希望它吐出向量时，我会吐出一个标量。这应该通过更改

create_I

:

中的 return 语句来解决

def create_I(f, a, b, N):
    def I(x):
        return trap_1D(lambda y: f(x, y), a, b, N)
    
    return np.vectorize(I)