标签: ode

这是我尝试过的一个循环,std::vector<double>并且使用普通的旧版本double*.

对于1000万个元素,矢量版本始终在该double*版本所占用的时间的80%左右运行; 几乎任何价值N,矢量明显更快.

偷看GCC STL源代码,我没有看到它std::vector正在做的事情比double*成语正在做的更好(即用普通的旧分配new[],operator[]取消引用偏移量). 这个问题也说明了这一点.

任何想法为什么矢量版本更快？

Compiler: GCC 4.6.1
Example compile line: g++ -Ofast -march=native -DNDEBUG \
                      -ftree-vectorizer-verbose=2 -o vector.bin \
                      vector.cpp -lrt
OS: CentOS 5
CPU: Opteron 8431
RAM: 128 GB

如果我使用icpc 11.1或在Xeon上运行,结果在质量上是相同的.此外,矢量化器转储表示只有std::vector向量化的构造函数中的填充操作.

矢量版本:

#include <vector>
#include <iostream>
#include <boost/lexical_cast.hpp>
#include "util.h"
#include "rkck_params.h"

using namespace std;

int main( int argc, char* argv[] )
{
    const size_t N = …

我有一个二阶微分方程，我想在 python 中解决它。问题是，对于其中一个变量，我没有初始条件，0而只有无穷大的值。谁能告诉我应该提供哪些参数scipy.integrate.odeint？能解决吗？

方程：

需要在时间方面找到 Theta。它的一阶导数在 处为零t=0。theta 未知，t=0但它在足够长的时间内变为零。其余的都是已知的。由于近似值I可以设置为零，因此消除了二阶导数，这将使问题更容易。

我正在尝试使用 deSolve 在 R 中解决一个简单的 ODE： dQ/dt = f(Q)*(P - E).整个过程是 Q 的时间序列。诀窍是 P 和 E 是数据本身的固定时间序列，因此 diff eq 仅在 Q 中有效。

用固定的时间步长迭代解决这个问题相对简单，但我试图找到一种在 R 中使用自适应时间步长的方法。因为 P 和 E 是离散的，连续的时间步长可能具有相同的 P 和 E 值，其中很好。我一直在玩 deSolve，但一直无法解决这个问题。理想情况下，我想使用标准的四阶 Runge-Kutta。

有任何想法吗？在 MATLAB 中做吗？

编辑可复制的例子。我希望能够使用可变时间步长 Runge-Kutta 4 方法进行此计算。我可以很容易地对固定时间步 rk4 进行编程，而不是自适应。

working <- structure(list(datetime = structure(c(1185915600, 1185919200, 
1185922800, 1185926400, 1185930000, 1185933600, 1185937200, 1185940800, 
1185944400, 1185948000, 1185951600), class = c("POSIXct", "POSIXt"
), tzone = "UTC"), p = c(0, 0, 0, 1.1, 0.5, 0.7, 0, 0, 1.3, 0, …

我正在尝试使用来自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]每个时间步长后如何传递回函数？

我正在尝试数值求解一个允许离散跳跃的 ODE。我正在使用 Euler 方法，并希望 Numba 的 jit 可以帮助我加快进程（现在脚本需要 300 秒才能运行，而我需要它运行 200 次）。

这是我简化的第一次尝试：

import numpy as np
from numba import jit

dt = 1e-5
T = 1
x0 = 1
noiter = int(T / dt)
res = np.zeros(noiter)

def fdot(x, t):
    return -x + t / (x + 1) ** 2

def solve_my_ODE(res, fdot, x0, T, dt):
    res[0] = x0
    noiter = int(T / dt)
    for i in range(noiter - 1):
        res[i + 1] = res[i] + dt * fdot(res[i], …

我有一个看起来像的系统

dn/dt=f(n,v)
dh/dt=g(h,v)

我想在流形上解决这个方程F(v,n,h)=0,一个非线性函数v.我尝试使用类似的东西v=fzero(@(x) F(x,n,h),0)在每个时间步骤解决歧管上的v的值.但这非常慢,而且我的系统是张弛振荡器,无法满足积分容差.如何在定义的流形上找到ODE的解F(v,n,h)=0？

我正在研究使用的常微分方程系统deSolve,并想知道是否有任何方法可以防止差分变量值低于零.我已经看过一些关于在矢量,数据框等中将负值设置为零的其他帖子,但由于这是一个生物模型(并且T细胞计数变为负值没有意义),需要阻止它开始发生,所以这些值不会扭曲结果,而不仅仅是替换最终输出中的底片.

我有一个包含 3 个边界条件的 3 个微分方程系统（从我相信的代码中可以很明显地看出）。我设法在 MATLAB 中用一个循环来解决它，以便在程序即将返回错误时一点一点地改变初始猜测而不终止程序。然而，在scipy's solve_bvp，我总是得到一些答案，虽然它是错误的。所以我一直在改变我的猜测（这一直在改变答案）并且我给出的数字与我从实际解决方案中得到的数字非常接近，但它仍然无法正常工作。代码是否存在其他问题，因此无法正常工作？我刚刚编辑了他们的文档代码。

import numpy as np
def fun(x, y):
    return np.vstack((3.769911184e12*np.exp(-19846/y[1])*(1-y[0]), 0.2056315191*(y[2]-y[1])+6.511664773e14*np.exp(-19846/y[1])*(1-y[0]), 1.696460033*(y[2]-y[1])))
def bc(ya, yb):
    return np.array([ya[0], ya[1]-673, yb[2]-200])
x = np.linspace(0, 1, 5)
#y = np.ones((3, x.size))
y = np.array([[1, 1, 1, 1, 1], [670, 670, 670, 670, 670], [670, 670, 670, 670, 670] ])
from scipy.integrate import solve_bvp
sol = solve_bvp(fun, bc, x, y)

实际解决方案如下图所示。

BVP 的 MATLAB 解决方案

我试图实现分步傅立叶方法来解决光学中的非线性薛定谔方程.它主要分别处理线性部分和非线性部分.它通过傅里叶变换和时域中的非线性部分来求解线性部分.

从书中复制以下代码:

alpha = 0
beta_2 = 1
gamma = 1

T = linspace(-5,5,2^13);
delta_T = T(2)-T(1);

L = max(size(A));
delta_omega = 1/L/delta_T*2*pi;
omega = (-L/2:1:L/2-1)*delta_omega;

A = 2*sech(T);
A_t = A;
step_num = 1000;
h = 0.5*pi/step_num;
results = zeros(L,step_num);

A_f = fftshift(fft(A_t));
for n=1:step_num
    A_f = A_f.*exp(-alpha*(h/2)-1i*beta_2/2*omega.^2*(h/2));
    A_t = ifft(A_f);
    A_t = A_t.*exp(1i*gamma*(abs(A_t).^2*h));
    A_f = fft(A_t);
    A_f = A_f.*exp(-alpha*(h/2)-1i*beta_2/2*omega.^2*(h/2));
    A_t = ifft(A_f);
    results(:,n) = abs(A_t);
end

A_t脉冲在哪里(要解决的功能).我不明白的是,它最初用于fftshift将零频率移动到中心,但后来又在循环中没有fftshift.我尝试添加fftshift到主循环,或在最开始删除它.两者都给出了错误的结果,为什么呢？在一般情况下,我应该什么时候使用fftshift和 …

这是我做的:

from sympy import *
x = symbols("x")
y = Function("y")
dsolve(diff(y(x),x) - y(x)**x)

我得到的答案(SymPy1.0)是:

Eq(y(x), (C1 - x*(x - 1))**(1/(-x + 1)))

但那是错的.双方Mathematica并Maple不能解决这个ODE.这里发生了什么事？