python - 使用黎曼和的中点规则计算 2 变量定积分的体积

Question

我正在尝试 sin^2(x)-cos^2(y)使用while和for循环来近似 2 变量定积分的体积。我经常更改代码，最近一次更改，它坏了。我对 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,0,1,0,2))
    
def Riemann(x0,xn,y0,yn,N):
    e = 1; 
    while e > 1e-3:
        x = np.linspace(0,1,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,1,0,2,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,1,0,2,10))
     22 
     23 

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

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

标签： pythonnumpyscipy

Answer 1

Your code has multiple issues, I'll address them so that you can improve. First off all, the formatting is pretty terrible. Put spaces after commas, separate things more with white-space, etc. You can see the difference in my code, and I'm by no means an expert at code formatting.

Second, your method is not doing what you think it's doing. Every time you iterate i, you create an entire array of values and you assign it to V, since x and y are both arrays. Neither x nor y are being updated here. The loop does the same thing every time, and V get re-assigned the same value every time. np.sum(V) never gets assigned anywhere, so the only thing getting updated at all in that loop is e. Of course, that bit is incorrect since you cannot subtract a vector from a scalar, since, as I wrote above, V is a vector.

Your function didn't use x0, y0, etc. for your bounds of integration, since your linspaces were hardcoded.

Now we come to the solution. There are two approaches to this problem. There's the "slow" pure Python way, where we just loop over our y's and x's and take function values multiplied by the base dx * dy. That version looks like this:

# This a naive version, using two for loops. It's very slow.
def Riemann(x0, xn, y0, yn, N):
    xs = np.linspace(x0, xn, N)
    ys = np.linspace(y0, yn, 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

Note I moved the multiplication outside to save some on performance.

The other way is to use numpy, that version looks like this:

def Riemann(x0, xn, y0, yn, N):
    points = itertools.product(np.linspace(x0, xn, N), np.linspace(y0, yn, 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))

Here we avoid the double for-loop. Note that you must include import itertools for this to work. Here we use the Cartesian product to create all points we wish to evaluate, and then give those points to your function f which is designed to work with numpy arrays. We get a vector back from f of all the function values at each point, and we simply just sum all the elements, just like we did in the for-loop. Then we can multiply by the common base dx * dy and return that.

The only thing I do not understand about your code is what you want e to do, and how it relates to N. I'm guessing it's some sort of error tolerance, but why you were trying to subtract the total volume so far (even if your code did nothing of the sort) from 1 I can't understand.