时间戳

首页 > 解决方案 > quad 和 trapz 积分返回非常不同的值

python - quad 和 trapz 积分返回非常不同的值

问题描述

我想使用quad 而不是trapz 来制作更正式的代码，但是quad 返回一个非常不同的值，并且还说“IntegrationWarning：已达到最大细分数（50）......”。这是代码：

def delta_t(m_lens,z_lens,y_impact):
return 4*const.G*const.c**-3*m_lens*const.M_sun*(1+z_lens)*(0.5*y_impact*np.sqrt(y_impact**2+4)+np.log((np.sqrt(y_impact**2+4)+y_impact)/(np.sqrt(y_impact**2+4)-y_impact)))

def y_max(r_max):
y_max = np.sqrt((1+r_max)/r_max**0.5-2)
return y_max

def y_min(m_lens,z_lens,dt_min):
y_min = fsolve(lambda y: delta_t(m_lens, z_lens, y) - dt_min*u.s, 0)[0]
return y_min

def cosmo_dependent_part(z_lens,z_source,f_dm):
cosmo_dependent_part = 3/2*f_dm*0.24/const.c*((cosmo.WMAP9.H0*1000*u.meter/u.kilometer)**2*cosmo.WMAP9.angular_diameter_distance_z1z2(0,z_lens)*cosmo.WMAP9.angular_diameter_distance_z1z2(z_lens,z_source))/(cosmo.WMAP9.H(z_lens)*1000*u.meter/u.kilometer*cosmo.WMAP9.angular_diameter_distance_z1z2(0,z_source))
return cosmo_dependent_part

第一个积分：

def optical_depth_fixed_source_integrand(z_source,z_lens,f_dm,m_lens,dt_min,r_max):
return cosmo_dependent_part(z_lens,z_source,f_dm)*(1+z_lens)**2*(y_max(r_max)**2-y_min(m_lens,z_lens,dt_min)**2)

def optical_depth_fixed_source(z_source,f_dm,m_lens,dt_min,r_max):
    z_lens_array = np.linspace(0,z_source,100)
    optical_depth_fixed_source_integrand_array=np.array([optical_depth_fixed_source_integrand(z_source,z_lens,f_dm,m_lens,dt_min,r_max) for z_lens in z_lens_array])
    optical_depth_fixed_source = integrate.trapz(optical_depth_fixed_source_integrand_array,z_lens_array)
    return optical_depth_fixed_source
print(optical_depth_fixed_source(1,1,30,1e-3,5))

输出：

0.017963036494289607

第二个积分：

def optical_depth_fixed_source0(z_source,f_dm,m_lens,dt_min,r_max):
optical_depth_fixed_source0 = integrate.quad(optical_depth_fixed_source_integrand,0,z_source,args=(z_source,f_dm,m_lens,dt_min,r_max))[0]
return optical_depth_fixed_source0
print(optical_depth_fixed_source0(1,1,30,1e-3,5))

第二个输出：

-5.17389415131513
C:\Users\USER\anaconda3\lib\site-packages\ipykernel_launcher.py:2: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.

什么是“专用集成商”，我该怎么做？

我尝试将四元积分拆分为较小的段，但没有帮助。

我还能尝试什么？

如何提高限​​额？

我无法最小化代码，所以我添加了所有代码

标签： pythonscipyintegration

解决方案

推荐阅读