标签: numerical-integration

我有一个358个数字的向量.我想对这个向量进行数值积分,但我不知道这个向量的功能.

我发现我们可以使用trapz或quad,但我真的不明白如何在没有函数的情况下进行集成.

在MATLAB中,ode45有一个参数调用NonNegative,它将解决方案约束为非负.他们甚至写了一篇关于这种方法是如何工作的文章,以及它不是什么愚蠢的东西,只要将y_i设置为0就会变成负数,因为这通常不起作用.

现在,MATLAB也有dde23解决延迟微分方程的问题,但NonNegative该积分器没有等效参数.

不幸的是,我任务是将延迟到使用解决了现有的ODE系统ode45与NonNegative启用.

有什么想法我应该继续吗？

编辑:

我不确定这是否有用,但......

我的系统的DDE部分看起来像:

 dx = 1/(1+y*z) - x;
 dy = (y*z)^2/(1+(y*z)^2) - y;
 dz = X - z;

其中X(第三个等式中的大写字母变量)是延迟版本x.然后,我通过添加几个方面的方程此DDE系统连结至现有的(以及更大)的ODE系统x和z,然后将合并的系统集成在一起.

在下面的Python中,我有五个函数包含在func我必须集成的数组中.代码调用使用f2py以下代码生成的外部Fortran模块:

import numpy as np
from numpy import cos, sin , exp
from trapzdv import trapzdv
def func(x):
    return np.array([x**2, x**3, cos(x), sin(x), exp(x)])

if __name__ == '__main__':
    xs = np.linspace(0.,20.,100)
    ans =  trapzdv(func,xs,5)
    print 'from Fortran:', ans
    print 'exact:', np.array([20**3/3., 20**4/4., sin(20.), -cos(20.), exp(20.)])

Fortran例程是:

      subroutine trapzdv(f,xs,nf,nxs,result)
          integer :: I
          double precision :: x1,x2
          integer, intent(in) :: nf, nxs
          double precision, dimension(nf) :: fx1,fx2
          double precision, intent(in), dimension(nxs) :: xs
          double precision, intent(out), dimension(nf) :: …

我正在尝试使用来自scipy的odeint解决二阶ODE。我遇到的问题是该函数隐式地与二阶项耦合，如简化的代码段所示（请忽略示例的假装物理学）：

import numpy as np
from scipy.integrate import odeint

def integral(y,t,F_l,mass):
    dydt = np.zeros_like(y)
    x, v = y
    F_r =  (((1-a)/3)**2 + (2*(1+a)/3)**2) * v # 'a' implicit 
    a  = (F_l - F_r)/mass

    dydt = [v, a]
return dydt


y0 = [0,5]
time = np.linspace(0.,10.,21)
F_lon = 100.
mass = 1000.

dydt = odeint(integral, y0, time, args=(F_lon,mass))

在这种情况下，我意识到可以代数求解隐式变量，但是在我的实际场景中，逻辑之间存在很多逻辑F_r，a并且代数运算的评估失败。

我相信可以使用MATLAB的ode15i函数来解决DAE ，但我尝试尽可能避免这种情况。

我的问题是-有办法解决python（最好是scipy）中的隐式ODE函数（DAE）吗？有没有更好的方法解决以上问题呢？

作为最后的选择，可以接受上a一个时间步长。dydt[1]每个时间步长后如何传递回函数？

我已经尝试过由纽曼（Newman）编写的计算物理学进行练习，并为自适应梯形规则编写了以下代码。当每张幻灯片的误差估计值大于允许值时，它将把该部分分成两半。我只是想知道我还能做些什么来使算法更有效。

xm=[]
def trap_adapt(f,a,b,epsilon=1.0e-8):    
    def step(x1,x2,f1,f2):
        xm = (x1+x2)/2.0
        fm = f(xm)
        h1 = x2-x1
        h2 = h1/2.0
        I1 = (f1+f2)*h1/2.0
        I2 = (f1+2*fm+f2)*h2/2.0
        error = abs((I2-I1)/3.0) # leading term in the error expression
        if error <= h2*delta:
            points.append(xm) # add the points to the list to check if it is really using more points for more rapid-varying regions
            return h2/3*(f1 + 4*fm + f2)
        else:
            return step(x1,xm,f1,fm)+step(xm,x2,fm,f2) 
    delta = epsilon/(b-a)
    fa, fb = f(a), f(b)  
    return step(a,b,fa,fb)

此外，我使用了一些简单的公式将其与Romberg积分进行比较，发现对于相同的精度，该自适应方法使用更多的点来计算积分。

仅仅是因为其固有的局限性吗？使用这种自适应算法代替Romberg方法有什么优势？有什么方法可以使其更快，更准确？

Johnathan Dummer在网上有一个常用的verlet-integration公式,叫做Time-Corrected Verlet.但是我已经阅读了几个论坛帖子,人们在某些条件下会得到奇怪或意想不到的结果.

Johnathan Dummer的公式:

x1 = x + (x – x0) * dt / dt0 + a * dt^2

还有一个stackoverflow答案,它说明了Dummer的时间校正公式被打破了,并且海报将他自己的推导作为正确的推导.

通过stackoverflow答案建议正确的公式

x1 = x + (x – x0) * dt / dt0 + a * dt * (dt + dt0) / 2

嗯,达默的公式真的坏了吗？如果是的话,海报的推导更好吗？

PS:Dummer x1 = x - x0 + a * dt^2在他的网站上使用verlet集成公式而不是正确的,这也很奇怪x1 = 2x - x0 + a * dt^2.

如何优化数值积分例程（与 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 …

我想通过使用 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

我想编写一个使用牛顿方法的程序:

估计这个积分的x:

其中X是总距离.

我有函数计算通过使用梯形方法进行数值积分到达一定距离所需的时间.不使用trapz.

function T = time_to_destination(x, route, n)
h=(x-0)/n;
dx = 0:h:x;
y = (1./(velocity(dx,route)));

Xk = dx(2:end)-dx(1:end-1); 
Yk = y(2:end)+y(1:end-1);

T = 0.5*sum(Xk.*Yk);
end

它通过一组数据点之间的三次样条插值的ppval获取其速度值.外推值不应该是可取的.

function [v] = velocity(x, route)
load(route);

if all(x >= distance_km(1))==1 & all(x <= distance_km(end))==1
    estimation = spline(distance_km, speed_kmph);
    v = ppval(estimation, x);
else
    error('Bad input, please choose a new value')
end
end

速度样条曲线的绘图如果您对此感兴趣,请评估:

dx= 1:0.1:65

现在我想编写一个函数,可以解决在给定时间之后行进的距离,使用没有fzero/fsolve的newton方法.但我不知道如何解决积分的上界.

根据微积分的基本定理,我认为积分的导数是积分内的函数,这是我试图重新创建的Time_to_destination /(1/velocity)我添加了我想要解决的常数到时间目的地所以它

(Time_to_destination - (输入时间))/(1 /速度)

不确定我是否正确行事.

编辑:重写我的代码,现在效果更好,但我对Newton Raphson的停止条件似乎没有收敛到零.我也试图从梯形积分(ET)实现错误,但不确定我是否应该打扰实现它.还可以在底部找到路径文件.

牛顿方法的停止条件和误差计算:

梯形误差估计:

Function x = distance(T, route) …

我有兴趣使用该scipy.integrate库为给定的初始点集成向量场（即找到流线）。由于矢量场是numpy.ndarray在计算网格上定义的对象，因此必须对网格点之间的值进行插值。有没有集成商处理这个问题？也就是说，如果我要尝试以下操作

import numpy as np
import scipy.integrate as sc
vx = np.random.randn(10,10)
vy = np.random.randn(10,10)
def f(x,t):
    return [vx[x[0],x[1]], vy[x[0],x[1]]] # which obviously does not work if x[i] is a float
p0 = (0.5,0.5)
dt = 0.1
t0 = 0
t1 = 1
t = np.arange(t0,t1+dt,dt)
sc.odeint(f,p0,t)

编辑 ：

我需要返回周围网格点的向量场的插值：

def f(x,t):
    im1 = int(np.floor(x[0]))
    ip1 = int(np.ceil(x[1]))
    jm1 = int(np.floor(x[0]))
    jp1 = int(np.ceil(x[1]))
    if (im1 == ip1) and (jm1 == jp1):
        return [vx[x[0],x[1]], vy[x[0],x[1]]]
    else: …