标签: numerical-integration

我想知道是否有人知道基于numpy/scipy的python包在数值上整合了一个复杂的数值函数在一个细分域(在我的特定情况下,一个由voronoi单元限制的2D域)？在过去,我使用了几个matlab文件交换包,但是如果可能的话,我希望保持在我当前的python工作流程中.matlab例程是

http://www.mathworks.com/matlabcentral/fileexchange/9435-n-dimensional-simplex-quadrature

用于正交和网格生成使用:

http://www.mathworks.com/matlabcentral/fileexchange/25555-mesh2d-automatic-mesh-generation

关于网格生成以及随后对网格的数值积分的任何建议都将受到赞赏.

我正在寻找一个超级数字正交函数.它应该具有以下三个属性:

• 自适应 - 它会自动调整采样点的密度以适应被积函数.这是绝对必要的,因为我的被积函数非常不均匀且计算成本很高.
• 矢量化 - 为了提高效率,它会在采样点列表上调用被积函数,而不是一次调用一个点.
• 能够处理向量值函数 - 向量值被积函数的所有组件都是同时计算的,无需额外成本,因此将所有组件分别集成是没有意义的.

另外,它应该是:

• 2D - 我想要计算的积分是平面区域上的双积分,我希望能够指定整个积分的整体(相对)容差,并让它适当地管理误差预算.

有人知道有这样一个函数的库吗？即使四个属性中的两个或三个也不会好.

我正在使用Python和SciPy,所以如果它已经与Python一起使用,那就是奖励.(但是我也可以编写胶水代码,让它在必要时调用我的被积函数.)

从我的SymPy输出我有下面显示的矩阵,我必须在2D中集成.目前我正在以元素方式进行,如下所示.此方法适用,但它变得太慢(用于sympy.mpmath.quad和scipy.integrate.dblquad)为我的真实案例(其中A,其职能是更大(见下面编辑):

from sympy import Matrix, sin, cos
import sympy
import scipy
sympy.var( 'x, t' )
A = Matrix([[(sin(2-0.1*x)*sin(t)*x+cos(2-0.1*x)*cos(t)*x)*cos(3-0.1*x)*cos(t)],
            [(cos(2-0.1*x)*sin(t)*x+sin(2-0.1*x)*cos(t)*x)*sin(3-0.1*x)*cos(t)],
            [(cos(2-0.1*x)*sin(t)*x+cos(2-0.1*x)*sin(t)*x)*sin(3-0.1*x)*sin(t)]])

# integration intervals
x1,x2,t1,t2 = (30, 75, 0, 2*scipy.pi)

# element-wise integration
from sympy.utilities import lambdify
from sympy.mpmath import quad
from scipy.integrate import dblquad
A_int1 = scipy.zeros( A.shape, dtype=float )
A_int2 = scipy.zeros( A.shape, dtype=float )
for (i,j), expr in scipy.ndenumerate(A):
    tmp = lambdify( (x,t), expr, 'math' )
    A_int1[i,j] = quad( tmp, (x1, x2), …

我想在python中使用多个cpus在数字上集成一个函数.我想做的事情如下:

from scipy.integrate import quad
import multiprocessing
def FanDDW(arguments):
  wtq,eigq_files,DDB_files,EIGR2D_files,FAN_files = arguments
  ...
  return tot_corr

# Numerical integration
def integration(frequency):
# Parallelize the work over cpus
  pool = multiprocessing.Pool(processes=nb_cpus)
  total = pool.map(FanDDW, zip(wtq,eigq_files,DDB_files,EIGR2D_files,FAN_files))
  FanDDW_corr = sum(total)
  return quad(FanDDW, -Inf, Inf, args=(zip(wtq,eigq_files,DDB_files,EIGR2D_files,FAN_files)))[0]

vec_functionint = vectorize(integration)
vec_functionint(3,arange(1.0,4.0,0.5))

"频率"也是一个全局变量(FanDDW(参数)外部).它是一个包含必须计算函数的位置的向量.我想四元组应该以聪明的方式选择频率.如何将它传递给FanDDW,知道它不应该在CPU之间分配,并且pool.map正是这样做的(这就是为什么我把它作为全局变量并且没有将它作为参数传递给定义的原因).

感谢您的任何帮助.

塞缪尔.

我有一个DataFrame力 - 位移数据.位移数组已设置为DataFrame索引,列是不同测试的各种力曲线.

如何计算完成的工作("曲线下面积")？

我看着numpy.trapz哪些似乎做了我需要的东西,但我认为我可以避免像这样循环遍历每一列:

import numpy as np
import pandas as pd 

forces = pd.read_csv(...)
work_done = {}

for col in forces.columns:
    work_done[col] = np.trapz(forces.loc[col], forces.index))

我希望DataFrame在曲线下创建一个新的区域而不是a dict,并且认为DataFrame.apply()或某些东西可能是合适的但不知道从哪里开始寻找.

简而言之:

1. 我可以避免循环吗？
2. 我DataFrame可以直接创建一项工作吗？

在此先感谢您的帮助.

早上好/下午好/晚上好,

我正在研究一个涉及四阶张量计算体积积分的Matlab脚本.设H(r,theta,phi)是我想要整合的函数.假设不能通过对r,theta和phi的简单操作来获得H.

我的问题是在Matlab中和我知道的任何其他代码一样:

All input functions must accept arrays and operate elementwise. The function FUN(X,Y,Z)
must accept arrays X, Y, Z of the same size and return an array of corresponding values.

这是来自Matlab 的integral3函数的实现.

如果我尝试这个简单的功能:

fun = @(X,Y,Z) X.*Y.*Z

完全没有问题,如果我将它集成在[0,1] x [0,1] x [0,1]上,我得到了正确的结果:

integral3(fun,0,1,0,1,0,1)

返回0.125,这是正确的.

问题在于,正如我所说,我不能用向量进行简单的计算来获得H,我不得不这样做或多或少做事:

function [result] = fun(x,y,z)
    sz = length(x);
    result  = zeros(1,sz);
    for i=1:sz
        result(i) = x(i)*y(i)*z(i);
    end
end

这个函数独立工作,返回与我之前介绍的另一个完全相同的结果.但是,当我尝试使用integral3时,我收到此错误:

Error using integral2Calc>integral2t/tensor (line 241)
Integrand output size does not match the input …

我有这个联合测试版发行版,我正在努力整合.

alpha1 = 400
beta1 = 26000
alpha2 = 410
beta2 = 26000

integrate(function(y) {
   sapply(y, function(y) {integrate(function(x) dbeta(x, shape1 = alpha1, 
     shape2 = beta1)*dbeta(y, shape1 = alpha2, shape2 = beta2), y,1)$value})
}, 0, 1)

这适用于较小的alpha和beta甚至达到这些值.但如果我在测试版中大得多,那么整合功能就会开始破坏.我一直把它与蒙特卡洛的整合进行比较

n_sim = 1000000   # number of simulations
y = rbeta(n_sim, shape1 = alpha2, shape2 = beta2)
C2 = (1-pbeta(y, shape1 = alpha1, shape2 = beta1))
mean(C2)

这对于较小的合理alpha和beta产生大致相同的答案.是的,我需要更大的测试版,我理解为什么我需要这么大的测试版.它来自大型数据集.有人知道为什么会这样吗？也许集成功能的内部工作原理？有工作吗？

我有以下问题。我有一个f使用 numpy 函数在 python 中定义的函数。该函数在正实数上平滑且可积。我想构造该函数的双重反导数（假设反导数在 0 处的值和斜率均为 0），以便我可以在任何小于 100 的正实数上对其进行评估。

fat的反导数的定义x：

integrate f(s) with s from 0 to x

fat的双重反导数的定义x：

integrate (integrate f(t) with t from 0 to s) with s from 0 to x

的实际形式f并不重要，因此为了方便起见，我将使用一个简单的形式。但请注意，即使我的示例具有已知的封闭形式，但我的实际函数却没有。

import numpy as np

f = lambda x: np.exp(-x)*x

我的解决方案是使用朴素数值积分将反导数构造为数组：

N = 10000
delta = 100/N

xs = np.linspace(0,100,N+1)

vs = f(xs)
avs = np.cumsum(vs)*delta
aavs = np.cumsum(avs)*delta

这当然有效，但它给了我数组而不是函数。但这不是一个大问题，因为我可以aavs使用样条曲线进行插值来获取函数并摆脱数组。

from scipy.interpolate import …

我希望能够使用numpys trapz函数计算跟随积分

numpy.trapz([-1, 1]) # returns 0

但我不想允许负面区域.有没有一种有效的方法来做到这一点,还是我必须寻找最小点并手动转换数组？

是否numpy.trapz(numpy.abs([-1, 1]))有意义？

如果我有一个随机数Z，它被定义为另外两个随机数 和 的总和，X那么Y的概率分布是和Z的概率分布的卷积。卷积基本上是分布函数乘积的积分。卷积中的积分通常没有解析解，因此必须使用基本求积算法来计算。在伪代码中：XY

prob_z(z) = integrate(lambda t: prob_x(t) * prob_y(z-t), -inf, inf)\n

举一个具体的例子，可以使用以下Python/Scipy代码计算Z正态分布变量X和对数正态分布变量的总和：Y

from scipy.integrate import quad\nfrom scipy.stats import norm, lognorm\nfrom scipy import log\n\nprob_x = lambda x: norm.pdf(x, 0, 1)  # N(mu=0, sigma=1)\nprob_y = lambda y: lognorm.pdf(y, 0.1, scale=10)  # LogN(mu=log(10), sigma=0.1)\ndef prob_z(z):\n    return quad(lambda t: prob_x(t)*prob_y(z-t), -inf, inf)\n

现在我想计算对数概率。天真的解决方案是简单地执行以下操作：

def log_prob_z(z):\n    return log(prob_z(z))\n

然而，这在数值上是不稳定的。在大约 39 个标准差之后，概率分布在数值上为 …