问题标签 [numerical-integration]

我正在尝试在 Python 中重现这个 Mathematica 程序：

它找到数值积分的根，并形成这些值的图。但是，我无法尝试运行。

当前尝试：

from scipy.integrate import quad from scipy import integration from scipy.optimize import fsolve import pylab as pl import numpy as np

我正在使用 R 的积分函数，但被积函数给我带来了一些麻烦（它有一个奇点）。考虑到这一点，

有没有比使用截止参数并手动调整它直到给出结果（不一定准确）更好的方法来计算 R 中的此类主值？

我想评估形式为 \int_{-\infty}^a \int_{-infty}^b sum_{i,j}^K a_i*a_j*x^i*y^j*exp(- x^2 - y^2 + x*y)dx dy

其中 a_i 和 a_j 是常数。由于积分是线性的，我可以互换求和和积分，但在这种情况下，我必须评估 K^2 积分，而且时间太长。在这种情况下，我会执行以下操作：

时间太长，所以我只想评估一个积分，但我不知道如何向量化 sum_{i,j}^K a_i*a_j*x^i*y^j*exp(-x^2 - y^2 + x*y)，即如何提供给integral2。| 将非常感谢任何帮助。

我想评估形式为 $$\int_{-\infty}^a \int_{-\infty}^b \sum_{i,j}^K a_ia_jx^iy^j\exp(-x^ 2 - y^2 + xy)dx dy $$

其中 $a_i$ 和 $a_j$ 是常量。由于积分是线性的，我可以互换求和和积分，但在这种情况下，我必须评估 $K^2$ 积分，而且花费的时间太长。在这种情况下，我会执行以下操作：

时间太长，所以我只想评估一个积分，但我不知道如何向量化 $\sum_{i,j}^K a_ia_jx^iy^j\exp(-x^2 - y^2 + xy)$，即如何提供给integral2。如果有任何帮助，我将不胜感激。

我一直在尝试更多地了解泛型编程，因为我认为我对此了解的不够多。所以我正在考虑如何实现我的一个程序的模板版本。我尝试使用的程序是一个数值积分器程序，用户在其中选择要使用的积分器（即 Euler、Runge Kutta 等），然后对他们选择的任何函数进行积分。我目前这样做的方法是拥有一个名为 Integrator 的抽象基类，以及几个实现集成方法的派生类。所以代码看起来像这样（还有很多事情要做，但这只是为了展示方法）。请注意，我为此使用 Qt 并声明了一个 Integrator *integrator; 在 MainWindow 类中。

使用 getIntegrator 本质上使用工厂方法

所以这个方法效果很好，积分器程序运行的很好。现在我知道模板函数是在编译时生成的，并且考虑到我正在使用运行时信息，你将如何使用模板来实现它？如果问题不清楚，我要问的是......在运行时给定用户选择，即使用哪个积分器，我如何使用模板方法调用正确的积分函数？

我有以下代码：

我想针对 和 的Og不同值和不同值来评​​估它xmax，ymax在这种特殊情况下是Og = [4 3;3 4], xmax = 3, ymax = 3。有什么办法可以在没有循环的情况下做到这一点？我想避免循环，因为我的速度有问题。我试过cellfun了，但它没有帮助。

从我的 SymPy 输出中，我得到了如下所示的矩阵，我必须将其集成为 2D。目前我正在按元素进行操作，如下所示。这种方法有效，但对于我的实际情况（其中和它的功能要大得多（见下面的编辑）来说，它变得太慢了（对于sympy.mpmath.quad和）：scipy.integrate.dblquadA

我正在考虑一次性完成，但我不确定这是否是要走的路：

编辑：真实案例已在此链接中提供。只需解压缩并运行shadmehri_2012.py（此示例的作者来自：Shadmehri et al. 2012）。我已经为可以执行以下操作的人开始了 50 的赏金：

• 使其比提出的问题合理地快
  
• 即使有许多术语m=15和n=15代码也能在不出现内存错误的情况下运行），我管理到m=732n=7位

当前时间可以总结如下（用 m=3 和 n=3 测量）。由此可见，数值积分是瓶颈。

构建试验函数 = 0%
评估微分方程 = 2%
羔羊化 k1 = 22%
积分 k1 = 74%
羔羊化和积分 k2 = 2%
提取特征值 = 0%

相关问题：关于lambdify

我想计算一个积分，它由两个函数确定： I(T) = ∫<sub>0 ^T i( f(t), g(t)) dt 其中f和g求解常微分方程，而i是已知。

显而易见的方法是推导出I的微分方程，然后将其与f和g一起求解（这可以完成，但在我的情况下在数值上很昂贵）。但是，在我的情况下，f求解具有初始条件f(0)和g的方程以及具有最终条件g(T)的方程。

目前我最好的猜测是使用标准 ODE 求解器在网格上求解f和g，然后使用标准方法进行数值积分，使用等距t坐标或某种正交规则（基本上是数值食谱中描述的任何内容） .

有没有人有更好的解决方案？也就是说，一种将特定类型的 ode 求解器及其精度考虑在内的方法。

我有一个函数 f(x) = 1/(x + a+ b*I*sign(x)) 我想计算积分

dx dy dz f(x) f(y) f(z) f(x+y+z) f(xy - z)

在整个 R^3 上（b>0 和 a,- b 是统一的）。这只是一个有代表性的例子——实际上我有 n<7 个变量和 2n-1 个 f() 实例，其中 n 个涉及 n 个积分变量，其中 n-1 个涉及积分变量的一些线性组合。在这个阶段，我只对相对误差为 1e-3 左右的粗略估计感兴趣。

我尝试了以下库：

• Steven Johnson 的 cubature 代码：hcubature 算法可以工作，但速度非常慢，即使 n=2 也需要数亿次被积函数。
• HintLib：我尝试了与 Genz-Malik 规则、体积例程、VEGAS 和 MISER 与 Mersenne twister RNG 的自适应集成。对于 n = 3，只有第一个似乎是可行的选择，但对于 n = 3 和 relerr = 1e-2，它再次需要数亿次被积函数评估，这并不令人鼓舞。

对于积分区域，我尝试了两种方法：在 [-200, 200]^n 上进行积分（即一个区域如此之大以至于它基本上捕获了大部分积分）和替换 x = sinh(t) 这似乎是标准把戏。

我在数值分析方面没有太多经验，但大概困难在于 sign() 术语的不连续性。对于 n=2 和 f(x)f(y)f(xy)，沿 x=0、y=0、x=y 存在不连续性。这些在原点周围产生了一个非常尖锐的峰值（在各个象限中具有不同的符号），并且在 x=0,y=0,x=y 处形成了一种“脊”，其中被积函数的绝对值很大并且符号变化为你越过他们。所以至少我知道哪些地区很重要。我在想也许我可以做蒙特卡洛，但以某种方式提前“告诉”算法将注意力集中在哪里。但我不太确定该怎么做。

如果您对如何以合理的计算能力评估积分或如何使我的蒙特卡洛“想法”发挥作用有任何建议，我将不胜感激。我已经坚持了一段时间，所以欢迎任何意见。提前致谢。

我尝试在 Google 上搜索库以在 CUDA 上进行数值集成，但找不到任何库。

1）我想问一下，是否有任何库可用于在 CUDA 上执行（函数的）集成？

2) 如果我在 CUDA 上编写自己的代码，例如实现 Romberg 集成，我应该如何进行？假设我有功能，比如说f(x)；我是否需要针对不同的区间计算此函数的积分，例如0.0 - 0.1, ..., 0.2 - 0.3, ..., 1.3 - 2.3？如何并行计算所有这些？

在我看来，策略是如果我必须执行，例如，1000集成，我生成1000线程，每个线程计算梯形以及错误估计。但是，如果我想计算其中一个积分区间的梯形以及其他积分，我不知道如何以编程方式处理这个问题。