关键词：
【中文标题】使用 Numba 求解 ODE【英文标题】：Solve ODEs with Numba 【发布时间】：2021-11-24 06:54:56 【问题描述】：我正在尝试使用 Numba 使我的 ODE 求解器更快，但以下代码会引发键入错误：
import numpy as np
import matplotlib.pyplot as plt
from numba import njit
@njit
def pend(t, y, b, c):
theta, omega = y
dydt = np.array([omega, -b*omega - c*np.sin(theta)])
return dydt
@njit
def rungeStep(f, t, y0, tau, params):
k1 = tau * f(t, y0, *params)
k2 = tau * f(t, y0 + k1 / 2, *params)
k3 = tau * f(t, y0 + k2 / 2, *params)
k4 = tau * f(t, y0 + k3, *params)
return (k1 + 2 * k2 + 2 * k3 + k4) / 6
@njit
def integrate(f, t0, y0, tEnd, h, params):
ys = y0.copy()
t = np.array(t0)
while t0 <= tEnd:
y0 += rungeStep(f, t0, y0[0], h, params)
t0 += h
ys = np.concatenate((ys, y0), axis=0)
t = np.append(t, t0)
return t, ys.T
args = (0.25, 5)
y0 = np.array([[np.pi - 0.1, 0.0]])
t, y = integrate(pend, 0, y0, 10, 1, args)
这会导致：
TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Cannot unify array(int64, 0d, C) and array(int64, 1d, C) for 't.2', defined at <ipython-input-56-38d2ea70b889> (6)
File "<ipython-input-56-38d2ea70b889>", line 6:
def inagrate(f, t0, y0, tEnd, h, params):
<source elided>
while t0 <= tEnd:
y0 += rungeStep(f, t0, y0[0], h, params)
^
During: typing of assignment at <ipython-input-56-38d2ea70b889> (6)
File "<ipython-input-56-38d2ea70b889>", line 6:
def inagrate(f, t0, y0, tEnd, h, params):
<source elided>
while t0 <= tEnd:
y0 += rungeStep(f, t0, y0[0], h, params)
^
没有njit-decorator 它可以正常工作。有人可以帮帮我吗？
【问题讨论】：
【参考方案1】：如果没有 JIT，应该会得到同样的错误，但广播规则可能足够灵活。你的y0 是一个二维数组，你传递给rungeStep（Heun-Kutta 方法）是一个一维数组，它的返回值也是如此。因此，您不能立即使用该结果更新y0，您必须更新y0[0]。但问题是，为什么首先要引入这种复杂性。
Numpy 数组被优化为固定大小的对象，附加到它们需要在每个实例中分配内存和复制，这比附加到列表要慢，其中重新分配是通过（增长的）保留完成的。 Numba-JIT 似乎在将 numpy 数组列表转换为 2d numpy 数组时存在问题。有效的方法是将数组手动降级为简单列表。经过反复试验，我运行了以下代码（仅重复有更改的部分）：
@njit
def rungeStep(f, t, y0, tau, params):
k1 = tau * f(t, y0, *params)
k2 = tau * f(t + tau/2, y0 + k1 / 2, *params)
k3 = tau * f(t + tau/2, y0 + k2 / 2, *params)
k4 = tau * f(t + tau, y0 + k3, *params)
return (k1 + 2 * k2 + 2 * k3 + k4) / 6
@njit
def integrate(f, t0, y0, tEnd, h, params):
ys = [list(y0)]
t = [t0]
while t0 <= tEnd:
y0 += rungeStep(f, t0, y0, h, params)
t0 += h
ys.append(list(y0))
t.append(t0)
return np.array(t), np.array(ys).T
args = (0.25, 5.0)
y0 = np.array([np.pi - 0.1, 0.0])
t, y = integrate(pend, 0, y0, 10, 1e-1, args)
请注意，您的 RK4 实现有一些遗漏，在这里并不重要，但对于依赖时间的 ODE 示例来说，这是一个错误源。
绘制结果给出合理的图形
【讨论】：
非常感谢！你真的救了我！您可能知道这是调整 ODE 求解器大小的最佳选择，还是有一些更好的（使用 numba 更快）选项（RK4 func 除外，因为我在实际问题中不需要自适应步骤）。我对 numba 真的不太了解。提前谢谢！ RK4 通常是“经典”简单的固定步长四阶方法。您的意思可能是 RK45，它是具有嵌入式 4 阶方法的 Dormand-Prince 5th，许多标准求解器使用它来实现具有自适应步长的方法。如果您关心速度，您可以使用 JITcode，它将 ODE 函数转换为 C 代码，并直接与 odeint 后面的 Fortran 方法接口，或 PyDStool，它是灵活而全面的 julia-lang diffeq 库的包装器。使用 SymPy 表达式和 SciPy 求解器求解一阶 ODE 系统
】使用SymPy表达式和SciPy求解器求解一阶ODE系统【英文标题】：SolvingafirstordersystemofODEsusingSymPyexpressionsandSciPysolver【发布时间】：2018-09-2018:06:35【问题描述】：我正在尝试学习如何将pythonSymPy与SciPy集成以数值求解常微分方程。但... 查看详情
使用 scipy 求解第一个 ODE 的运动方程
】使用scipy求解第一个ODE的运动方程【英文标题】：SolvemotionequationsforfirstODEusingscipy【发布时间】：2021-12-3014:00:35【问题描述】：我想使用scipysolve_ivp函数求解运动一阶ODE方程。我可以看到我做错了什么，因为这应该是一个椭圆，... 查看详情
如何使用 scipy.integrate.odeint 求解具有时间相关变量的 ODE 系统
】如何使用scipy.integrate.odeint求解具有时间相关变量的ODE系统【英文标题】：HowtosolveasystemofODEswithscipy.integrate.odeintwithatime-dependentvariable【发布时间】：2017-08-1608:30:51【问题描述】：我正在使用scipy食谱中的ZombieApocalypseexample来学习... 查看详情
为啥 MATLAB 在使用 ode 求解器时会更改矩阵维度？
】为啥MATLAB在使用ode求解器时会更改矩阵维度？【英文标题】：WhydoesMATLABchangethematrixdimensionswhenusinganodesolver?为什么MATLAB在使用ode求解器时会更改矩阵维度？【发布时间】：2012-11-0119:42:03【问题描述】：我已经为一些运行良好的... 查看详情
deSolve：无法理解如何使用根函数提前停止 ode 求解器
】deSolve：无法理解如何使用根函数提前停止ode求解器【英文标题】：deSolve:Can\'tunderstandhowtoearlystopodesolverwithrootfunctions【发布时间】：2022-01-1909:55:24【问题描述】：我对如何在满足特定条件时停止求解器感到困惑。我准备了一个... 查看详情
在 numba / NumbaLSODA 中使用 eval()？
】在numba/NumbaLSODA中使用eval()？【英文标题】：useeval()innumba/NumbaLSODA?【发布时间】：2021-11-2402:44:51【问题描述】：我正在尝试使用NumbaLSODA来优化我的模型，这实质上是在解决初始值ODE问题（可能非常僵硬）。我的原始模型基于sci... 查看详情
matlab使用高阶求解器解决天体力学问题
使用ode78和ode89解决天体力学问题，该问题需要ODE求解器的每一步都具有高精度才能成功积分。和ode45都ode113无法使用默认容错来解决问题。即使使用更严格的错误阈值，ode89由于它在每个步骤中使用的Runge-Kutta公式的高精度，在... 查看详情
通过 Python 并行求解具有大量初始条件的 ODE
...【发布时间】：2020-01-0212:56:42【问题描述】：我正在尝试使用scipy.integrate.ode或scipy.integrate.odeint为一大组（超过一千个）初始条件求解ODE系统，但是执行循环非常慢，scipy似乎没有提 查看详情
在基类和派生类中分离 ODE 和 ODE 求解器
...有两个类A和B。B继承自A。B中有一个成员函数（add）需要使用A中的成员函数运行。classA;typedefint(A::*operation)(intc);t 查看详情
求解隐式 ODE（微分代数方程 DAE）
...【发布时间】：2014-06-2800:27:27【问题描述】：我正在尝试使用scipy中的odeint解决二阶ODE。我遇到的问题是函数隐式耦合到二阶项，如简化的sn-p所示（请忽略示例的假装物理）：importnumpyasnpfromscipy.integrateimporto 查看详情
matlab求解具有多个初始条件的ode系统
...具有多组初始条件的常微分方程组的方法。这些技术是：使用for-loop执行多个模拟，每组初始条件一个。这种技术使用简单，但不能为大型系统提供最佳性能。向量化ODE函数以同时求解所有初始条件集的方程组。这种技术对于大... 查看详情
关于 PDE，我可以在每个步骤中使用 ODE 求解器吗？
】关于PDE，我可以在每个步骤中使用ODE求解器吗？【英文标题】：InmatterofPDEs,canIuseODEsolversforeachstep?【发布时间】：2021-12-3006:49:53【问题描述】：我有这样的方程式：wl[n]=w[n]+Δx*v[n]ϕl[n]=ϕ[n]+Δx*ρ[n]ρl[n]=ρ[n]-Δt*fρ(ρ,v,w,n)vl[n]=v[n]-... 查看详情
matlab求解刚性ode
什么是刚性ODE？对于某些ODE问题，与积分区间相比，求解器采用的步长被强制降低到不合理的小水平，即使在解曲线平滑的区域也是如此。这些步长可能非常小，以至于遍历一个短时间间隔可能需要数百万次评估。这可能导致求... 查看详情
matlab求解具有强状态相关质量矩阵的ode
使用移动网格技术[1]求解Burgers方程。该问题包括一个质量矩阵，并且指定选项以解决质量矩阵的强状态依赖性和稀疏性，从而使求解过程更加高效。问题设置Burgers方程是由PDE给出的对流扩散方程∂u∂t=ϵ∂2u∂x2−∂∂x 查看详情
如何用ode45求解matlab中的耦合微分方程(代码片段)
...。答案这很容易。首先编写一个函数来实现微分方程，并使用与函数名对应的文件名保存：functiondy=my_ode(t,y)dy(1)=y(1)*(0.3/y(1)^3+0.)^(1/2);%ady(2)=1/dy(1);%tau然后在MATLAB中，使用您的函数调用ode45求解器[t,y]=ode45(@my_ode,[010],[1;0]);这是结果... 查看详情
使用最小平方的非线性 ODE 优化
】使用最小平方的非线性ODE优化【英文标题】：nonlinearODEsoptimizationwithleastsq【发布时间】：2021-11-3006:28:55【问题描述】：[更新]我正在研究非线性ODE系统优化并将其拟合到实验数据。我有一个由5个模型ODE组成的系统，必须通过17... 查看详情
matlab求解非刚性ode
非刚性vanderPol方程vanderPol方程是二阶ODE在哪里 查看详情
时变输入ode求解(代码片段)
我正在尝试使用ode45解决一组ODE方程。我的一些参数已经是时间的函数，但我不断收到错误。functionododx(1,1)=(vl+vr)/2*cos(x(3));dx(2,1)=(vl+vr)/2*sin(x(3));dx(3,1)=obz其中obz，vr和vl都是例如：obz=sin(t),t=0:dt:tf;我使用以下语法：[t,x1]=ode45(@(t,x)od... 查看详情