在 python 中向量化二阶 ode 求解

计算科学 Python 一体化 矢量化
2021-12-23 20:41:39

我正在尝试编写一个 python 程序,通过数值积分二阶常微分方程来模拟大量粒子的运动。我首先将 ODE 拆分为两个耦合的一阶 ODE,并使用与该问题scipy.integrate.solve_ivp的答案中描述的类似方法进行求解但是,我的问题是我想为大量粒子解决这个系统,每个粒子都有不同的初始条件。天真地我可以用一个循环来做到这一点,但我确信 numpy 和 scipy 必须有一种方法可以更快地矢量化这个操作。for

我查看了文档scipy.integrate.solve_ivp它谈到了矢量化,但不是我想要的方式。我希望对于具有两个耦合的一阶 ODE 和 n 个粒子的系统,您可以将初始条件作为具有大小的数组输入,(2,n)但事实并非如此。

有没有办法解决多个初始条件而不求助于慢 pythonfor循环?

作为参考,我想解决的 ODE 系统看起来像

dx/dt = v
dv/dt = F(x,v)

在数组中具有初始条件,例如

initialConditions = [[x0,v0],[x1,v1],...,[xN,vN]]
1个回答

有没有办法解决多个初始条件而不求助于慢 pythonfor循环?

不幸的是,我认为您的问题的解决方案可能是并行化,而不是矢量化。虽然 Python比 C 或 C++ 慢,但所for讨论的循环实际上并不是计算的一部分,由于初始条件不同,它在数值上是一个不同的问题。此外,计算中涉及的数据太多,无论如何缓存行上都没有多余的空间。

就矢量化而言,我敢打赌,在这方面没有什么可以做的来改进代码。集成算法在 Fortran 中实现,由于其严格的别名规则,它甚至比 C 和 C++ 快很多倍,允许更积极的优化

就 Python 代码本身而言,尤其是您的循环,您甚至可能很难通过用 C 编写它来显着提高执行速度,因为这部分代码的影响可以忽略不计与集成期间实际完成的繁重工作相比,这甚至没有考虑将 Python 转换为 C 需要多长时间。

好处是 AVX-512 指令集已经出现很久了的,甚至可能是一些可靠的 GPU 来启动。它可能会产生足够的影响,你永远不知道;不幸的是,优化问题的定义总是非常具体,所以大多数问题的答案通常只是“取决于情况”。抱歉,我无法提供更多帮助。

通常,这里的标准建议是编写尝试运行一些试验,分析执行情况,确定瓶颈在哪里,然后将所有精力集中在这上面,因为即使是几个百分点的改进也可能意味着几十个(甚至几百个)数小时的宝贵计算时间。不幸的是,scipy 代码可能并没有变得更好,所以在这种情况下,我可能会尝试以某种方式减少要进行分析测试的初始条件的绝对数量。例如,如果您正在使用 Navier-Stokes 方程,我正在考虑诸如平均之类的东西。