使用黎曼求和的中点法则计算 2 变量定积分的体积

如何解决使用黎曼求和的中点法则计算 2 变量定积分的体积

我正在尝试使用 sin^2(x)-cos^2(y) 和 while 循环来近似计算 2 个可变定积分 for 的体积。我经常更改代码，最近一次更改，它坏了。我对 python 很陌生，所以我仍在弄清楚如何正确使用数组。

这就是我现在所拥有的（编辑：通过 alani 的评论，我设法修复了错误，但现在我在运行代码时没有收到答案）

import numpy as np
import scipy.integrate


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

print(scipy.integrate.dblquad(f,1,2))
    
def riemann(x0,xn,y0,yn,N):
    e = 1; 
    while e > 1e-3:
        x = np.linspace(0,N)
        y = np.linspace(0,2,N)
        dx = (x0-xn)/N
        dy = (y0-yn)/N
        for i in range(N):
            V = (dx*dy)*(f(x,y))
            np.sum(V)
            e = abs(1-V)
print(riemann(0,1000))

运行此代码时，我收到：

(-0.2654480895858587,9.090239973208559e-15)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-9-c654507b2f73> in <module>
     19             np.sum(V)
     20             e = abs(1-V)
---> 21 print(riemann(0,10))
     22 
     23 

<ipython-input-9-c654507b2f73> in riemann(x0,N)
     10 def riemann(x0,N):
     11     e = 1;
---> 12     while e > 1e-3:
     13         x = np.linspace(0,N)
     14         y = np.linspace(0,N)

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

解决方法

您的代码有多个问题，我会解决这些问题以便您改进。首先，格式非常糟糕。逗号后面加空格，用空格分隔更多的东西，等等。你可以看到我的代码的不同之处，我绝不是代码格式方面的专家。

其次，您的方法没有按照您认为的那样做。每次迭代 i 时，都会创建一个完整的值数组并将其分配给 V，因为 x 和 y 都是数组。 x 和 y 均未在此处更新。循环每次都做同样的事情，V 每次都被重新分配相同的值。 np.sum(V) 永远不会被分配到任何地方，因此在该循环中完全唯一得到更新的是 e。当然，这个位是不正确的，因为你不能从标量中减去一个向量，因为正如我上面写的，V 是一个向量。

您的函数没有使用 x0、y0 等作为您的集成边界，因为您的 linspace 是硬编码的。

现在我们来解决这个问题。有两种方法可以解决这个问题。有一种“慢”的纯 Python 方式，我们只是循环遍历 y 和 x 并取函数值乘以基数 dx * dy。那个版本看起来像这样：

# This a naive version,using two for loops. It's very slow.
def Riemann(x0,xn,y0,yn,N):
    xs = np.linspace(x0,N)
    ys = np.linspace(y0,N)
    dx = (x0 - xn) / N
    dy = (y0 - yn) / N
    V = 0
    for y in ys:
        for x in xs:
            V += f(x,y)
    return dx * dy * V

请注意，我将乘法移到外面以节省一些性能。

另一种方式是使用 numpy，那个版本看起来像这样：

def Riemann(x0,N):
    points = itertools.product(np.linspace(x0,N),np.linspace(y0,N))
    points = np.array(list(points))
    xs = points[:,0]
    ys = points[:,1]
    dx = (x0 - xn) / N
    dy = (y0 - yn) / N
    return dx * dy * np.sum(f(xs,ys))

这里我们避免了双重 for 循环。请注意，您必须包含 import itertools 才能使其正常工作。在这里，我们使用笛卡尔积来创建我们希望评估的所有点，然后将这些点提供给您的函数 f，该函数旨在处理 numpy 数组。我们从每个点的所有函数值的 f 返回一个向量，我们只是简单地对所有元素求和，就像我们在 for 循环中所做的那样。然后我们可以乘以公共基数 dx * dy 并返回它。

关于您的代码，我唯一不明白的是您希望 e 做什么，以及它与 N 的关系。我猜这是某种容错能力，但我不明白为什么到目前为止你试图从 1 中减去总体积（即使你的代码没有做任何类似的事情）。