问题描述
我正在尝试使用 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 中减去总体积(即使你的代码没有做任何类似的事情)。