标签: numerical-integration

我正在尝试在Python 2.7.2中实现梯形规则.我写了以下函数:

def trapezoidal(f, a, b, n):
    h = float(b - a) / n
    s = 0.0
    s += h * f(a)
    for i in range(1, n):
        s += 2.0 * h * f(a + i*h)
    s += h * f(b)
    return s

但是,f(lambda x:x**2,5,10,100)返回583.333(它应该返回291.667),所以很明显我的脚本有问题.我不能发现它.

我花了几个小时研究数值积分和速度/位置估计，但我真的找不到一个可以被我的大脑理解或适合我的情况的答案。

我有一个 IMU（惯性测量单元），它有一个陀螺仪、一个加速度计和一个磁力计。所有这些传感器都处于融合状态，这意味着例如使用陀螺仪我能够补偿加速度计读数中的重力，而磁力计补偿漂移。换句话说，我可以使用这样的设置获得纯加速度读数。

现在，我正在尝试根据加速度准确估计位置，您可能知道这需要双重积分，并且有多种方法可以做到这一点。但我不知道哪一个最适合这里。有人可以分享一些有关此的信息吗？另外，如果您能在不使用任何复杂的数学公式/符号的情况下向我解释它，我将不胜感激，我不是数学家，这是我在寻找信息时遇到的问题之一。

谢谢

朱莉娅绝对是新人，有一个问题要问你。

我正在通过移植我的一些 Mathematica 和 Python 代码（主要是物理学中的科学计算等）来自学一些 Julia，然后看看是什么。到现在为止，事情还算顺利。而且很快。到目前为止。

现在，我正在模拟一个基本的锁定放大器，它本质上接受一个可能非常复杂的时间相关信号Uin(t)，并产生一个输出，Uout(t)在某个参考频率上锁相fref（也就是说，它突出显示的分量Uin(t)，与参考正弦波有一定的相位关系）。描述无关紧要，重要的是它基本上是通过计算积分来实现的（为了清楚起见，我实际上在这里省略了阶段）：

因此，我在 Mathematica 和 Julia 中开始并测试了这一点：我定义了一个模型Uin(t)，传递了一些参数值，然后Uout(t)在时间t = 0为fref.

• Julia：我使用QuadGK包进行数值积分。

  T = 0.1
  f = 100.
  Uin(t) = sin(2pi * f * t) + sin(2pi * 2f *t)
  Uout(t, fref)::Float64 = quadgk(s -> Uin(s) * sin(2pi * fref * s), t-T, t, rtol=1E-3)[1]/T
  frng = 80.:1.:120.
  print(@time broadcast(Uout, 0., frng))
  

  Run Code Online (Sandbox Code Playgroud)

• 数学

  T = …

  Run Code Online (Sandbox Code Playgroud)

我需要整合一个在其表达式中具有阶乘的函数。但是，如果您尝试计算阶乘，当 n > 170 时，R 将返回 Inf。

我发现很多包可以让你计算非常大的数字，但是，它们总是从我无法集成的类中返回一个对象。积分的最终结果总是一个小数。

这是我的代码：

integrand <- function(n, i, x) {
    (factorial(n) / (factorial(i - 1) * factorial(n - i))) *
        x^(i - 1) * (1 - x)^(n - i)
}

forder <- function(Fx, x, n, i, ...) {
    lower <- sapply(x - 1, Fx, ...)
    upper <- sapply(x, Fx, ...)
     integrate(integrand,
               lower = lower, upper = upper, n = n, i = i,
               stop.on.error = FALSE)$value
}
forder <- Vectorize(forder, "x")

##------------------------------------------------------------------------------
## Some example …

我有以下 Matlab ODE 代码：

[t,y,~,~,ie] = ode23tb(@(t,y) RHSODE(t,y),[0,t_end], [i0;v0],options);

我希望 ODE 求解器还可以给出结果 z，它是 y 和 dy/dt 的函数，这样 z = f(y,dy/dt)。

有谁知道如何将这样的 z 添加到求解器的输出中？

当集成普通的pdf时,我从scipy integrate.quad函数得到了一些奇怪的输出.这是我正在尝试使用的代码:

inpdf = lambda c: norm.pdf(50, loc=c, scale = 1)
result = integrate.quad(inpdf, -np.inf, np.inf)
print result

哪个返回(3.281718223506531e-99,0.0)这显然是错误的.当我将x值更改为0时,在下面的代码段中,我得到了适当的输出:

inpdf = lambda c: norm.pdf(0, loc=c, scale = 1)
result = integrate.quad(inpdf, -np.inf, np.inf)
print result

返回(0.9999999999999998,1.0178191320905743e-08)非常接近1,这应该是什么.使用R的积分函数完全相同的事情,所以这可能只是正交算法的一个已知属性？有谁知道为什么会这样？

我在R中的集成功能遇到了一些问题.我试图绘制积分的声音,但似乎我做得不好.

t <- seq(0, 0.04, 0.0001)
vi <- function(x) {5 * sin(2 * pi * 50 * x)}
vo <- function(x) {integrate(vi, lower=0, upper=x)$value}

test_vect = Vectorize(vo, vectorize.args='x')
plot(t, vo(t))  # should be a cosine wave
plot(t, vi(t))  # sine wave

vo应该是一个正弦波,但使用test_vect给我错误的情节和vo直接使用给出错误'x'和'y'长度不同.请问有人可以帮我解决这个问题吗？

我很清楚如何以scipy.integrate.ode.integrate(t)最简单的形式使用该函数,但API读取它还需要两个可选参数,即step和relax.当前文档没有关于这些参数的信息,也没有在示例中使用它们.我想知道,他们做了什么以及他们有用的情况是什么？

采取以下Mathematica代码：

f[x_] := Exp[-x];
c = 0.9;
g[x_] := c*x^(c - 1)*Exp[-x^c];
SetPrecision[Integrate[f[x]*Log[f[x]/g[x]], {x, 0.001, \[Infinity]}],20]

Mathematica可以毫无问题地进行计算并给出答案0.010089328699390866240。我希望能够执行类似的积分，但是没有Mathematica的副本。例如，仅凭天真地在scipy中实现它，就很难使用标准正交库，因为f（x）和g（x）任意接近0。这是使用标准正交的Python示例，由于需要无限的精度而失败：：

from scipy.integrate import quad
import numpy as np

def f(x):
    return sum([ps[idx]*lambdas[idx]*np.exp(- lambdas[idx] * x) for idx in range(len(ps))])

def g(x):
    return scipy.stats.weibull_min.pdf(x, c=c)

c = 0.9
ps = [1]
lambdas = [1]
eps = 0.001  # weibull_min is only defined for x > 0
print(quad(lambda x: f(x) * np.log(f(x) / g(x)), eps, np.inf)) # Output 

应该大于0

如何在代码中像Mathematica一样执行这种不正确的积分？我不介意使用哪种免费语言/图书馆。

我正在实现一个 RKF4(5) 积分器，我无法确定我的代码是否正常工作并且我不明白本地截断错误，或者我的代码是否无法正常工作。

对于代码块的大小，我深表歉意，但在这种情况下，最小可重现示例相当大。

import numpy as np

def RKF45(state, derivative_function, h):
    """
    Calculate the next state with the 4th-order calculation, along with a 5th-order error
    check.

    Inputs:
    state: the current value of the function, float
    derivative_function: A function which takes a state (given as a float)
                         and returns the derivative of a function at that point
    h: step size, float
    """
    k1 = h * derivative_function(state)
    k2 = h * derivative_function(state + (k1 / 4))
    k3 = h * …