用scipy积分一个复杂函数

kx7yvsdv  于 2023-05-07  发布在  其他
关注(0)|答案(1)|浏览(138)

我想对一个函数求数值积分。在这个函数中,我有两种类型的球坐标(r,th)和(rp,thp)。我想在(rp,thp)上积分,并尝试使用if循环为(r,th)给予。

import numpy as np
from scipy.integrate import quad

wavelength = 1e-3  # m
R = 5  # m radius
aperture_diameter = 5e-1  # m
k = 2 * np.pi / wavelength  # the wavenumber

# Define the number of points in the ring
num_points = 100
r_values = np.linspace(0, aperture_diameter / 2, num_points)
th_values = np.linspace(0, 2 * np.pi, num_points)

def f(rp, thp):
    integrand = np.zeros_like(rp)
    for i, rp_val in enumerate(rp):
        for j, th_val in enumerate(thp):
            integrand[i] = rp_val * np.exp(-1j * k / (2 * R) * (r_values[i]**2 + rp_val**2 - 2 * r_values[i]**2 * rp_val**2 * (np.cos(th_values[j]) * np.cos(th_val) + np.sin(th_values[j]) * np.sin(th_val))))
    return integrand

# Define the integration limits
rp_limits = [0, int(aperture_diameter / 2)]
thp_limits = [0, int(2 * np.pi)]

# Perform numerical integration
result, _ = quad(lambda rp, thp: f(rp, thp), *rp_limits, *thp_limits)

print(result)

但我得到了错误

TypeError: 'float' object is not iterable

请帮帮我,我很困惑。

jyztefdp

jyztefdp1#

你真的读过,并试图理解错误吗?还是跑来求我们帮忙
以下是完整的错误消息:

In [44]: result, _ = quad(lambda rp, thp: f(rp, thp), *rp_limits, *thp_limits)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[44], line 1
----> 1 result, _ = quad(lambda rp, thp: f(rp, thp), *rp_limits, *thp_limits)

File ~\miniconda3\lib\site-packages\scipy\integrate\_quadpack_py.py:411, in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    408 flip, a, b = b < a, min(a, b), max(a, b)
    410 if weight is None:
--> 411     retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
    412                    points)
    413 else:
    414     if points is not None:

File ~\miniconda3\lib\site-packages\scipy\integrate\_quadpack_py.py:523, in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    521 if points is None:
    522     if infbounds == 0:
--> 523         return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    524     else:
    525         return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

Cell In[44], line 1, in <lambda>(rp, thp)
----> 1 result, _ = quad(lambda rp, thp: f(rp, thp), *rp_limits, *thp_limits)

Cell In[43], line 15, in f(rp, thp)
     13 def f(rp, thp):
     14     integrand = np.zeros_like(rp)
---> 15     for i, rp_val in enumerate(rp):
     16         for j, th_val in enumerate(thp):
     17             integrand[i] = rp_val * np.exp(-1j * k / (2 * R) * (r_values[i]**2 + rp_val**2 - 2 * r_values[i]**2 * rp_val**2 * (np.cos(th_values[j]) * np.cos(th_val) + np.sin(th_values[j]) * np.sin(th_val))))

TypeError: 'float' object is not iterable

所以是enumerate(rp)出了问题。quad传递给frp是一个浮点数。你在期待什么?
在之前的追踪中

quad(func, a, b, args,...

所以func就是你的f(lambda不会改变任何东西)。ab来自*rp_limits*thp_limits被接收为args(或者可能只是thp_limits[0])。
在任何情况下,你似乎写f和限制,好像你希望做一个2变量的积分。正如quad签名(以及其文档的其余部分)应该清楚说明的那样,quad只集成了一个变量。
使用quad或其他类似的scipy函数时,请确保输入符合文档中的规格。它们应该清楚地告诉您func需要什么参数和其他各种参数。quad具有完全控制权;你就得顺应它的行为。

相关问题