符号演算与Python的集成
我试着用数值计算一个双积分。问题是，(我想)我需要符号积分和数值积分的混合。积分看起来是这样的：
我不能使用numpy.integrate，因为它不只是一个双积分，因为它的幂(1/a)在中间。我无法得到最里面的积分(然后提高到幂)的数字，因为它最终是一个依赖于x的函数，然后我需要对它进行积分。
我尝试使用符号演算，使用嵌套的sym.integrate，如下所示
sym.integrate((sym.integrate(sym.exp(-(w**2)/(2*sigmaw)-alpha*((x-w)**2)/(2*sigma)),(w,-sym.oo, sym.oo)))**(1/alpha),(x,-sym.oo, sym.oo))
然而，它只是吐回表达式本身，而没有数字。我想我需要得到一个内部积分的符号表达式，作为数值积分的函数。
有可能吗？如果不是在python中，使用另一种语言，如R？
对这类东西有经验吗？
回答 3
Stack Overflow用户
发布于 2022-01-19 05:36:18
我使用极大值( Maxima，https://maxima.sourceforge.io)，因为OP似乎在说所使用的确切系统并不太重要。
被积体只是高斯凸起的乘积，所以它在实线上的积分并不太难。Maxima没有世界上最强的集成商，但无论如何，它似乎可以处理这个问题。
首先假设所有参数都是正的；如果没有指定，Maxima将在计算期间请求符号。
(%i2) assume (alpha > 0, sigmaw > 0, sigma > 0);
(%o2) [alpha > 0, sigmaw > 0, sigma > 0]定义内部被积体。
(%i3) I: exp(-(w**2)/(2*sigmaw)-alpha*((x-w)**2)/(2*sigma));
2 2
alpha (x - w) w
(- --------------) - --------
2 sigma 2 sigmaw
(%o3) %e计算内部积分。
(%i4) I1: integrate (I, w, minf, inf);
(%o4) (sqrt(2) sqrt(%pi) sqrt(sigma) sqrt(sigmaw)
2
alpha x
- ------------------------
2 alpha sigmaw + 2 sigma
%e )/sqrt(alpha sigmaw + sigma)漂亮的打印机(ASCII艺术)在这里很难读懂，也许这种一维表示更有意义。grind产生一维显示.
(%i5) grind(%);
(sqrt(2)*sqrt(%pi)*sqrt(sigma)*sqrt(sigmaw)
*%e^-((alpha*x^2)/(2*alpha*sigmaw+2*sigma)))
/sqrt(alpha*sigmaw+sigma)$
(%o5) done定义外部被积体。
(%i7) I2: I1^(1/alpha);
1 1 1 1
------- ------- ------- -------
2 alpha 2 alpha 2 alpha 2 alpha
(%o7) (2 %pi sigma sigmaw
2
x 1
- ------------------------ -------
2 alpha sigmaw + 2 sigma 2 alpha
%e )/(alpha sigmaw + sigma)计算外积分。最后的结果在这里名为foo。
(%i9) foo: integrate (I2, x, minf, inf);
1 1 1 1
------- + 1/2 ------- ------- -------
2 alpha 2 alpha 2 alpha 2 alpha
(%o9) (%pi 2 sigma sigmaw
1
-------
2 alpha
sqrt(2 alpha sigmaw + 2 sigma))/(alpha sigmaw + sigma)计算参数的特定值的外部积分。
(%i10) ev (foo, alpha = 3, sigma = 3/7, sigmaw = 7/4);
1/6 1/6 1/6 1/3 2/3
2 3 7 159 %pi
(%o10) ----------------------------
sqrt(14)
(%i11) float(%);
(%o11) 5.790416728790489计算一个数值近似。注quad_qagi适用于无限间隔。
(%i12) ev (quad_qagi (lambda([x], quad_qagi (I, w, minf, inf)[1]^(1/alpha)), x, minf, inf),
alpha = 3, sigma = 3/7, sigmaw = 7/4);
(%o12) [5.790416728790598, 7.216782674725913E-9, 270, 0]看起来这支持符号结果。
(%i13) first(%) - %o11;
(%o13) 1.092459456231154E-13外部积分，在一维显示中，可能对复制到另一个程序很有用：
(%i14) grind(foo);
(%pi^(1/(2*alpha)+1/2)*2^(1/(2*alpha))*sigma^(1/(2*alpha))
*sigmaw^(1/(2*alpha))
*sqrt(2*alpha*sigmaw+2*sigma))
/(alpha*sigmaw+sigma)^(1/(2*alpha))$
(%o14) done如果可能的话，我强烈建议尝试获得一个象征性的结果；数值积分通常是很棘手的。在给出的例子中，如果你只能做内部积分，而不能做外部积分，那仍然是一个很大的胜利。你可以把内部积分的符号解插入到外部积分的数值近似中。
Stack Overflow用户
发布于 2022-01-18 17:19:59
数值积分是通过在参数的某些值上采样被积体来工作的。特别是牛顿-科特斯公式的样本均匀，而不同风格的高斯积分样本不规则.
所以在你的例子中，积分器将需要对x的各种值的内部积分进行评估，从而在x上进行积分，这意味着每次在w上与已知的x进行数值积分。
请注意，由于您的域是无界的，您必须使用变量的更改来使其有限。
如果内部积分有一个解析表达式，你当然可以使用它，并在x上进行数值积分。
Stack Overflow用户
发布于 2022-01-18 17:29:38
这不能回答你的问题，但它肯定会对你有帮助，因为其他人已经指出了其他有用的工具。
对于手头的集成，您实际上不需要进行象征性的集成。
数值积分仅仅是在定义的有限网格上求和，而w上的积分只是在w轴上求和，与x相同。
主要的问题是如何选择积分网格，因为它不可能是无限的，对于高斯人，我会说至少是他们的西格玛的10倍，因为你可以得到尽可能低的误差，至于网格间距，我会使它尽可能小，你可以等待它运行。因此，对于上面的集成，这将是等价的，确保您没有增加网格步骤，直到您有多少内存需要图片，否则您的pc将挂起。
import numpy as np
# define constants
sigmaw = 0.1
sigma = 0.1
alpha = 0.2
# define grid
max_w = 2
min_w = -max_w
min_x = -3
max_x = -min_x
steps_w = 2000 # don't increase this too much or you'll run out of memory
steps_x = 1000 # don't increase this too much or you'll run out of memory
dw = (max_w - min_w) / steps_w
dx = (max_x - min_x) / steps_x
x_vec = np.linspace(min_x, max_x, steps_x)
w_vec = np.linspace(min_w, max_w, steps_w)
x, w = np.meshgrid(x_vec, w_vec, sparse=True)
# do integration
inner_term = np.exp(-(w ** 2) / (2 * sigmaw) - alpha * ((x - w) ** 2) / (2 * sigma))
inner_integral = np.sum(inner_term, axis=0) * dw
del inner_term # to free some memory
inner_integral_powered = inner_integral ** (1 / alpha)
del inner_integral # to free some memory
outer_integral = np.sum(inner_integral_powered) * dx
print(outer_integral)https://stackoverflow.com/questions/70758801
复制相似问题
腾讯云开发者
Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有
深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 
粤公网安备44030502008569号
腾讯云计算（北京）有限责任公司 京ICP证150476号 | 京ICP备11018762号