标签: numerical-integration

如何对无限范围内的一维积分进行数值积分(采用什么数值方法,以及使用什么技巧),其中被积函数中的一个或多个函数是1d量子谐振子波函数.其中我想在谐振子基础上计算某些函数的矩阵元素:

披_Ñ(X)= N _ñ ħ _Ñ(X)EXP(-x ² /2)
其中H _Ñ(x)是厄米多项式

V _m,n =\int _ { - infinity} ^ {infinity} phi _m(x)V(x)phi _n(x)dx

同样在存在具有不同宽度的量子谐波波函数的情况下.

问题是波函数phi _n(x)具有振荡行为,这对于大n是一个问题,并且来自GSL(GNU Scientific Library)的自适应Gauss-Kronrod积分算法花费很长时间来计算,并且具有大的误差.

我想在一个维度上执行数值积分,其中被积函数是向量值的.integrate()只允许标量积分,因此我需要多次调用它.该cubature软件包看起来非常适合,但它似乎对1D积分表现不佳.考虑以下示例(标量值积分和1D积分),

library(cubature)
integrand <- function(x, a=0.01) exp(-x^2/a^2)*cos(x)
Nmax <- 1e3
tolerance <- 1e-4

# using cubature's adaptIntegrate
time1 <- system.time(replicate(1e3, {
  a <<- adaptIntegrate(integrand, -1, 1, tol=tolerance, fDim=1, maxEval=Nmax)
}) )

# using integrate
time2 <- system.time(replicate(1e3, {
  b <<- integrate(integrand, -1, 1, rel.tol=tolerance, subdivisions=Nmax)
}) )

time1
user  system elapsed 
  2.398   0.004   2.403 
time2
user  system elapsed 
  0.204   0.004   0.208 

a$integral
> [1] 0.0177241
b$value
> [1] 0.0177241

a$functionEvaluations
> [1] …

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

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

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

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

我们如何numpy.polynomial.legendre.leggauss在除了[-1, 1]？之外的时间间隔内使用NumPy包？

以下示例与scipy.integrate.quad间隔中的Gauss-Legendre方法进行比较[-1, 1].

import numpy as np
from scipy import integrate

# Define function and interval
a = -1.
b =  1.
f = lambda x: np.cos(x)

# Gauss-Legendre (default interval is [-1, 1])
deg = 6
x, w = np.polynomial.legendre.leggauss(deg)
gauss = sum(w * f(x))

# For comparison
quad, quad_err = integrate.quad(f, a, b)

print 'The QUADPACK solution: {0:.12} with error: {1:.12}'.format(quad, quad_err)
print 'Gauss-Legendre solution: {0:.12}'.format(gauss)
print 'Difference between …

在不使用质量矩阵的情况下，ode 求解器ode45可以求解 y'=f(t,y)。

但是对于涉及“质量”矩阵的问题，在 ode 求解器中有一个质量矩阵选项，M(t,y)y'=f(t,y)。

“质量”矩阵究竟是什么？这个术语是否来自质量弹簧系统振荡的质量？我在文档中找不到关于此的示例代码。此外，似乎我可以在 y'=f(t,y) 等式中的 f(t,y) 中对有关 t 和 y 的信息进行编码。在什么情况/示例中会出现 M(t,y)y'=f(t,y) 需要 M(t,y) 的地方？

如果我有一个随机数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 个标准差之后，概率分布在数值上为 …

我有一些代码使用 scipy.integration.cumtrapz 来计算采样信号的反导数。我想使用辛普森规则而不是梯形。然而 scipy.integration.simps 似乎没有累积的对应物......我错过了什么吗？有没有一种简单的方法可以与“scipy.integration.simps”进行累积集成？

如何优化数值积分例程（与 C 相比）？

目前已经做了什么：

1. 我用未装箱的向量替换了列表（显而易见）。
2. 我应用了“Read World Haskell”一书中描述的分析技术http://book.realworldhaskell.org/read/profiling-and-optimization.html。我内联了一些琐碎的函数，并在各处插入了很多刘海。这带来了大约 10 倍的加速。
3. 我重构了代码（即提取iterator函数）。这带来了 3 倍的加速。
4. 我尝试用 Floats 替换多态签名，如这个问题 Optimizing numeric array Performance in Haskell 的答案。这几乎提高了 2 倍的速度。
5. 我这样编译 cabal exec ghc -- Simul.hs -O2 -fforce-recomp -fllvm -Wall
6. 更新按照 cchalmers 的建议，type Sample = (F, F)被替换为 data Sample = Sample {-# UNPACK #-} !F {-# UNPACK #-} !F

现在的性能几乎和C代码一样好。我们可以做得更好吗？

{-# LANGUAGE BangPatterns #-}

module Main
  where

import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as UM …

是否可以在不使用 cubature 包的情况下在 R 中进行三重集成？

基于这篇文章中的答案

InnerFunc = function(x) { x + 0.805 }
InnerIntegral = function(y) { sapply(y, 
    function(z) { integrate(InnerFunc, 15, z)$value }) }
integrate(InnerIntegral , 15, 50)
16826.4 with absolute error < 1.9e-10

例如，要编码这个三重积分：

我试过

InnerMostFunc = function(v) { v + y^2 }
InnerMostIntegral = function(w) { sapply(w, 
   function(x) { integrate(InnerMostFunc, 1, 2)$value }) }
InnerIntegral = function(y) { sapply(y, 
   function(z){integrate(InnerMostIntegral, 1, 2)$value }) }
integrate(InnerIntegral, 0, 1)

我想通过使用 R 中的集成函数将 exp(-x) 从 0 积分到 100000。但是我发现答案是 2.061453e-45，几乎是 0（零）。真正的答案是 1-exp(-100000)，几乎为 1。如何使用 R 中的集成函数进行积分以接近正确的解决方案？

以下是使用的R代码

ab<-function(x) { return(exp(-x)) }
integrate(ab,0,100000)$value

输出是

 2.061453e-45