python scipy integrate_python – 使用scipy.integrate.quad来整合复数
把它分成实部和虚部是什么问题? scipy.integrate.quad需要集成的函数返回浮动(也称为实数)用于其使用的算法。import scipyfrom scipy.integrate import quaddef complex_quadrature(func, a, b, **kwargs):def real_func(x):return scipy.real(func(x))def i
把它分成实部和虚部是什么问题? scipy.integrate.quad需要集成的函数返回浮动(也称为实数)用于其使用的算法。
import scipy
from scipy.integrate import quad
def complex_quadrature(func, a, b, **kwargs):
def real_func(x):
return scipy.real(func(x))
def imag_func(x):
return scipy.imag(func(x))
real_integral = quad(real_func, a, b, **kwargs)
imag_integral = quad(imag_func, a, b, **kwargs)
return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
例如。,
>>> complex_quadrature(lambda x: (scipy.exp(1j*x)), 0,scipy.pi/2)
((0.99999999999999989+0.99999999999999989j),
(1.1102230246251564e-14,),
(1.1102230246251564e-14,))
这是你期望舍入误差 – exp(ix)从0的整数,pi / 2是(1 / i)(e ^ i pi / 2-e ^ 0)= -i(i-1)= 1 〜(0.99999999999999989 0.99999999999999989j)。
并且为了记录,万一它不是所有人都清楚,整合是一个线性的功能,意味着∫{f(x)kg(x)} dx =∫f(x)dx k∫g(x)dx (其中k是相对于x的常数)。或者对于我们的具体情况∫z(x)dx =∫Re z(x)dx i∫Im z(x)dx as z(x)= Re z(x)i Im z(x)
如果您尝试在复杂平面中的路径(除了沿着实际轴线)或复杂平面中的区域进行集成,则需要一个更复杂的算法。
注意:Scipy.integrate不会直接处理复杂的集成。为什么?在FORTRAN QUADPACK库中,特别是在qagse.f中,它明显要求功能/变量在实现“基于每个子间隔内的21点高斯 – 克罗诺德正交”的全局自适应正交之前,由Peter Wynn的epsilon加速,算法。”所以除非你想尝试修改底层的FORTRAN以使其处理复杂的数字,将它编译成一个新的库,否则你不会让它工作。
如果您真的想要在完全一个集成中使用复数的Gauss-Kronrod方法,请查看wikipedias page并直接实现如下(使用15-pt,7-pt规则)。注意,我记住函数重复常见的变量的公共调用(假设函数调用很慢,就像函数非常复杂一样)。也只有7点和15点规则,因为我不喜欢自己计算节点/权重,那些是维基百科上列出的,而是在测试用例(〜1e-14)中得到合理的错误,
import scipy
from scipy import array
def quad_routine(func, a, b, x_list, w_list):
c_1 = (b-a)/2.0
c_2 = (b+a)/2.0
eval_points = map(lambda x: c_1*x+c_2, x_list)
func_evals = map(func, eval_points)
return c_1 * sum(array(func_evals) * array(w_list))
def quad_gauss_7(func, a, b):
x_gauss = [-0.949107912342759, -0.741531185599394, -0.405845151377397, 0, 0.405845151377397, 0.741531185599394, 0.949107912342759]
w_gauss = array([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830050505119, 0.279705391489277,0.129484966168870])
return quad_routine(func,a,b,x_gauss, w_gauss)
def quad_kronrod_15(func, a, b):
x_kr = [-0.991455371120813,-0.949107912342759, -0.864864423359769, -0.741531185599394, -0.586087235467691,-0.405845151377397, -0.207784955007898, 0.0, 0.207784955007898,0.405845151377397, 0.586087235467691, 0.741531185599394, 0.864864423359769, 0.949107912342759, 0.991455371120813]
w_kr = [0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525, 0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728, 0.204432940075298, 0.190350578064785, 0.169004726639267, 0.140653259715525, 0.104790010322250, 0.063092092629979, 0.022935322010529]
return quad_routine(func,a,b,x_kr, w_kr)
class Memoize(object):
def __init__(self, func):
self.func = func
self.eval_points = {}
def __call__(self, *args):
if args not in self.eval_points:
self.eval_points[args] = self.func(*args)
return self.eval_points[args]
def quad(func,a,b):
''' Output is the 15 point estimate; and the estimated error '''
func = Memoize(func) # Memoize function to skip repeated function calls.
g7 = quad_gauss_7(func,a,b)
k15 = quad_kronrod_15(func,a,b)
# I don't have much faith in this error estimate taken from wikipedia
# without incorporating how it should scale with changing limits
return [k15, (200*scipy.absolute(g7-k15))**1.5]
测试用例:
>>> quad(lambda x: scipy.exp(1j*x), 0,scipy.pi/2.0)
[(0.99999999999999711+0.99999999999999689j), 9.6120083407040365e-19]
我不相信错误估计 – 当从[-1到1]整合时,我从wiki获取了推荐的错误估计值,并且这些值对我来说似乎不合理。例如,与真相相比的错误是〜5e-15不〜1e-19。我确定如果有人咨询了num食谱,你可以得到一个更准确的估计。 (可能需要通过(a-b)/ 2来重复一些权力或类似的东西)。
回想一下,python版本比只是调用scipy的基于QUADPACK的集成两次不太准确。 (如果需要,你可以改进它)。

GitCode 天启AI是一款由 GitCode 团队打造的智能助手,基于先进的LLM(大语言模型)与多智能体 Agent 技术构建,致力于为用户提供高效、智能、多模态的创作与开发支持。它不仅支持自然语言对话,还具备处理文件、生成 PPT、撰写分析报告、开发 Web 应用等多项能力,真正做到“一句话,让 Al帮你完成复杂任务”。
更多推荐
所有评论(0)