问题标签 [odeint]

我尝试在 Mac OS X 10.9.2 g++ 5.1 上的 boost_1_55_0 中运行 [odeint 复杂状态类型示例代码。

下面的代码是解决 Stuart-Landau 振荡器的网站上的副本

上面的示例代码无法编译。9个错误生成并以：

什么是错误？

我正在研究使用最优控制方法（轨迹生成和优化）的无限维优化算法。我想对其进行轨迹优化的系统是非线性的，形式为 $\dot{x}(t) = f(x(t), u(t), t)$。换句话说，我有随时间变化的非线性输入。

是否可以使用 boost::odeint 解决这样的 ODE？我在文档中没有找到任何提示，但我可能只是没有看到它。

我正在尝试编写一个脚本来解决ODE in python和绘制结果。我正在使用 scipy.integrate.odeint 来完成这项任务。我遵循了一个简单的教程并修改了我的代码以使用我想要解决的 ODE

这段代码对于像 -x 这样的简单函数运行完美，但问题是它不适用于我在那里使用的函数。我可以从这种方法中得到一个情节，但无论我使用，我都会得到相同的情节

y=(1/(p*V*C)) * (A*e*b*np.power(T, 4) - (A*e*b*np.power(x, 4)) - ((1.32*A)/(D**0.25))*(np.power((x-T), 1.25)))

或者

y=(1/(p*V*C)) * (A*e*b*pow(T, 4) - (A*e*b*pow(x, 4)) - ((1.32*A)/(D**0.25))*(pow((x-T), 1.25)))

或者

y=(1/(p*V*C)) * (A*e*b*np.power(T, 4) - (A*e*b*np.power(x, 4))

我想即使这也给了我同样的情节

y=(1/(p*V*C)) * ((-A*e*b*np.power(x, 4))

此外，我的 e 值应该在 0 和 1 之间，但我必须使用巨大的数字来获得看起来像指数衰减的东西。我基本上是在尝试这样做：sfu.ca/~rca10/rants/convection.pdf。有问题的 ODE 位于第 4 页的顶部。链接中的人可以使用正常的发射率，但我不能，即使正在用相同的值解决相同的 ODE（我正在做完全相同的实验室）

我到底做错了什么？

所以我试图解决三个 ODE 的系统，并开发了以下代码来使用 ODEint 解决它们。但是当我运行时，ODEint 在调用我的方程组函数时遇到了问题。

在我的输出中，我得到：

更改 y0 和 t 的顺序仍然会导致相同的错误。我不明白这个错误，也没有引导我解决这个问题。任何帮助将不胜感激。

I'm trying to solve a set of coupled differential equations using scipy.integrate.odeint. However when I try to run the program I get the following error:

TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe' odepack.error: Result from function call is not a proper array of floats.

Here is the code I use:

Could it be because the equations from afleiden are calculated using Sympy and thus probably sypmpy expressions? If so, is there anything i can do about it? I tried using lambdify but that didn't work out.

我的第一个模块有以下三个模块：

我的第二个：

我的第三个：

现在，当我尝试从第三个主模块运行程序时出现以下错误：in func a_x = -(220**2)*x/N - (G * MBH * x)/NT TypeError: can't将序列乘以“元组”类型的非整数

所以我认为程序认为 a_x 是一个元组，我在这里要做的是用 odeint 的起始值 init_val 计算 a_x、a_y 和 a_z。

这段代码在没有 a_x 中的 (G*MBH*x)/NT 部分之前可以完美运行所以我认为我们必须弄清楚为什么程序认为 a_x 是一个序列。

所以基本上我使用 odeint 来解决我拥有的 ode 并且我希望将它用作另一个方程中的变量（最终实际上是我的解的方程）。

所以在这里我解决了第一个 ode（我的代码实际上是在解决一个 ODE 矩阵，但同样的问题仍然适用）

我的 f() 函数定义如下：

所以在所有这些好东西之后，我得到了一个数组或数字，点，我可以绘制以显示解决方案等。但我的问题是，我需要使用这个解决方案（为了简单起见，假设它是一个而不是解决方案矩阵）插入我实际上试图找到的功能：

其中 bernoulli() 只是一个返回该订单的伯努利数的函数，而 adjointOp() 只是给定订单的伴随运算符（交换子关系）。

我首先在 Mathematica 中做到了这一点，其中 Omeg 的解决方案是 InterpolationFunction，但这只是作为变量工作并相乘等。但是，我不知道如何在 python 中处理这种情况。

更新：

我已经在其他系统上运行了这个示例。在 Intel i7-3630QM、Intel HD4000 和 Radeon HD 7630M 上，所有结果都相同。对于 i7-4700MQ / 4800MQ，当使用 OpenCL 或 64 位 gcc 与 32 位 gcc 时，CPU 的结果是不同的。这是默认情况下使用 SSE 的 64 位 gcc 和 OpenCl 以及使用 387 数学的 32 位 gcc 的结果，当设置 -mfpmath=387 时，至少 64 位 gcc 会产生相同的结果。所以我必须阅读更多内容并尝试使用 x86 浮点。谢谢大家的答案。

我已经针对不同 OpenCL 设备上的十个系统运行了“编程 CUDA 和 OpenCL：使用现代 C++ 库的案例研究”中的 Lorenz 系统示例，并且得到了不同的结果：

1. Quadro K1100M (NVIDIA CUDA)

   r => xyz
   0.100000 => -0.000000000000000000000000 5.64444444444444444444444444444444444444444444444452
   11.1888890
   =
   > 5.212534 5.212534 5.212530
   

   22.277779 => 3.178553 2.579687 17.946903
   27.822224 => 5.008720 7.753564 16.377680
   33.366669 => -13.381100 -15.252210 36.107887
   38.911114 => 4.256534 6.813675 ​​23.838787 44.455555
   => -11.083726 0.691549 53.632290 50.000000
   => -8.624105 -15.728293 32.516193
   

2. Intel(R) HD Graphics 4600 (Intel(R) OpenCL)

   R => xyz
   0.100000 => -0.000000 -0.000000 0.000000
   5.644444 => -3.519253 -3.519250 4.644451
   11.188890 => 5.212531 5.212538 10.188890
   16.733334 => 6.477320 6.477326 15.733339
   

   22.277779 => 7.246771 7.398651 20.735369
   27.822224 => -6.295782 -10.615027 14.646572
   33.366669 => -4.132523 -7.773201 14.292910
   38.911114 => 14.183139 19.582197 37.943520 44.455555
   => -3.129006 7.564254 45.736408 50.000000
   => -9.146419 -17.006729 32.976696
   

3. Intel(R) Core(TM) i7-4800MQ CPU @ 2.70GHz (Intel(R) OpenCL)

   R => xyz
   0.100000 => -0.000000 -0.000000 0.000000
   5.644444 => -3.519254 -3.519251 4.644453
   11.188890 => 5.212513 5.212507 10.188900
   16.733334 => 6.477303 6.477296 15.733332
   

   22.277779 => -8.295195 -8.198518 22.271002
   27.822224 => -4.329878 -4.022876 22.573458
   33.366669 => 9.702943 3.997370 38.659538
   38.911114 => 16.105495 14.401397 48.537579 44.455555
   => -12.551083 -9.239071 49.378693 50.000000
   => 7.377638 3.447747 47.542763
   

如您所见，三个设备在 R=16.733334 的值上达成一致，然后开始出现分歧。

我在没有 VexCL 的情况下使用 odeint 运行了相同的区域，并获得了接近 CPU 运行时 OpenCL 的结果的结果：

香草味：

示例代码可以在这里找到：https ://github.com/ddemidov/gpgpu_with_modern_cpp/blob/master/src/lorenz_ensemble/vexcl_lorenz_ensemble.cpp

我不确定我在这里看到了什么？由于 CPU 结果彼此如此接近，这看起来像是 GPU 的问题，但由于我是 OpenCL 新手，我需要一些指示如何找到其根本原因。

我已经实现了以下方程组的解决方案

首先在 Mathematica 中，然后在 Python 中。

我在 Mathematica 中的代码：

从中我得到以下情节：Plot1（如果它给出 403 Forbidden 消息，请在 url 字段内按 enter）

后来我把同样的代码写进了python：

这是我得到的情节： Plot2（如果它给出 403 Forbidden 消息，请在 url 字段中按 enter）

如果在较小的字段中再次绘制 Mathematica 解：

他将得到与 python 解决方案类似的结果。只有轴'会错位。

为什么剧情差别这么大？我究竟做错了什么？

我怀疑我对模型的 python 实现是错误的，尤其是在计算 f1 的地方。或者，在这种情况下， plot() 函数对于绘制参数方程可能根本不方便。

谢谢。

ps：很抱歉没有拍到文字里面的图片，让你的生活变得艰难；我还没有足够的声望。

我正在使用 c++ boost/odeint 来集成化学反应网络。一切正常，除非我消除了一些反应。然后有时会在输出中出现-nan。这是该输出的一小部分。第一列是时间，其他列是物种的浓度。

第二个物种只依赖于自己，第七个物种是恒定的。

所以我的问题是：odeint 什么时候给出 -nan 作为输出？也许某处有一个 0/​​0 的划分，但所有方程只包含产品，例如

dxdt[5] = - c[2] * x[5] + c[3] * x[9] * x[3];

非常感谢您的任何建议。