尝试使用solve_ivp的矢量化选项，奇怪的是它会抛出一个错误，即y0必须是一维的。微量元素：
from scipy.integrate import solve_ivp
import numpy as np
import math
def f(t, y):
theta = math.pi/4
ham = np.array([[1,0],[1,np.exp(-1j*theta*t)]])
return-1j * np.dot(ham,y)
def main():
y0 = np.eye(2,dtype= np.complex128)
t0 = 0
tmax = 10**(-6)
sol=solve_ivp( lambda t,y :f(t,y),(t0,tmax),y0,method='RK45',vectorized=True)
print(sol.y)
if __name__ == '__main__':
main()
调用签名是 fun(t, y)。这里 t 是一个标量，ndarray y 有两个选项：它可以具有形状 (n,)；那么 fun 必须返回形状为 (n,) 的 array_like。或者它可以具有形状（n，k）；那么 fun 必须返回一个形状为 (n, k) 的 array_like，即每一列对应于 y 中的单个列。两个选项之间的选择由向量化参数决定（见下文）。矢量化实现允许通过有限差分更快地逼近雅可比行列式（刚性求解器所需）。
错误 ：
ValueError：
y0必须是一维的。
Python 3.6.8
scipy。版本 “1.2.1”
I am facing a problem I do not manage to solve. I would like to use nlme or nlmODE to perform a non linear regression with random effect using as a model the solution of a second order differential equation with fixed coefficients (a damped oscillator).
I manage to use nlme with simple models, but it seems that the use of deSolve to generate the solution of the differential equation causes a problem. Below an example, and the problems I …
我正在学习Haskell,并且正在尝试尽可能快地编写代码.在本练习中,我正在为一个简单的一维物理系统编写一个Euler积分器.
-O3.编译的.它运行在1.166秒.-O3.编译.它运行时间为21.3秒.-O3 -fllvm,它运行4.022秒.那么,我是否遗漏了一些优化我的Haskell代码的东西？
PS.:我使用了以下参数:1e-8 5.
C代码:
#include <stdio.h>
double p, v, a, t;
double func(double t) {
return t * t;
}
void euler(double dt) {
double nt = t + dt;
double na = func(nt);
double nv = v + na * dt;
double np = p + nv * dt;
p = np;
v = nv; …我正在尝试用来ode45解决ODE系统:
[X,Y]= ode45(@sys,[0, T],y0);
哪里,
function dy = sys(t,y)
dy(1) = f_1(y)
dy(2) = f_2(y)
dy(3) = f_3(y)
end
问题是该函数ode45需要y0是初始值[y_1(0), y_2(0), y_3(0)],而在我的系统中,我只有值[y_2(0), y_3(0), y_3(T)]可用.
从数学上讲,这组初始/终端条件应该足以确定系统,但是有什么方法可以用ode45MATLAB中的任何其他函数来处理它？
谢谢!
Julia的新手,试图测试ODE求解器的速度.我在教程中使用了Lorenz方程
using DifferentialEquations
using Plots
function lorenz(t,u,du)
du[1] = 10.0*(u[2]-u[1])
du[2] = u[1]*(28.0-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;1.0;1.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob,reltol=1e-8,abstol=1e-8,saveat=collect(0:0.01:100))
加载软件包的开始时间约为25秒,代码在Jupyter笔记本电脑的Windows 10四核笔记本电脑上运行了7秒.我知道Julia需要先预编译包,这就是加载时间这么长的原因吗？我发现25秒难以忍受.此外,当我使用不同的初始值再次运行求解器时,运行时间(~1s)要少得多,为什么会这样？这是典型的速度吗？
我目前正在尝试使用SciPy的integrate.ode包来解决一对耦合的一阶ODE:比如Lotka-Volterra捕食者 - 食饵方程.但是,这意味着在集成循环期间,我必须更新我在每次迭代时发送给方法的参数,并且只是跟踪前一个值并调用set_f_params()每次迭代似乎没有做到这一点.
hprev = Ho
pprev = Po
yh = np.zeros(0)
yp = np.zeros(0)
while dh.successful() and dp.successful() and dp.t < endtime and dh.t < endtime:
hparams = [alpha, beta, pprev]
pparams = [delta, gamma, hprev]
dh.set_f_params(hparams)
dp.set_f_params(pparams)
dh.integrate(dh.t + stepsize)
dp.integrate(dp.t + stepsize)
yh = np.append(yh, dh.y)
yp = np.append(yp, dp.y)
hprev = dh.y
pprev = dp.y
我在每次迭代时设置的值set_f_params似乎没有传播到回调方法,这并不是非常令人惊讶,因为网上没有任何示例似乎涉及"实时"变量传递给回调,但这是我能想到将这些值放入回调方法的唯一方法.
有没有人对如何使用SciPy数字整合这些ODE有任何建议？
使用numba.jit加快了右手端计算odeint从scipy.integrate正常工作：
from scipy.integrate import ode, odeint
from numba import jit
@jit
def rhs(t, X):
return 1
X = odeint(rhs, 0, np.linspace(0, 1, 11))
但是这样使用integrate.ode：
solver = ode(rhs)
solver.set_initial_value(0, 0)
while solver.successful() and solver.t < 1:
solver.integrate(solver.t + 0.1)
装饰器产生以下错误@jit：
capi_return is NULL
Call-back cb_f_in_dvode__user__routines failed.
Traceback (most recent call last):
File "sandbox/numba_cubic.py", line 15, in <module>
solver.integrate(solver.t + 0.1)
File "/home/pgermann/Software/anaconda3/lib/python3.4/site-packages/scipy/integrate/_ode.py", line 393, in integrate
self.f_params, self.jac_params)
File "/home/pgermann/Software/anaconda3/lib/python3.4/site-packages/scipy/integrate/_ode.py", line …我目前有一个具有时间依赖常数的颂歌系统.例如
def fun(u, t, a, b, c):
x = u[0]
y = u[1]
z = u[2]
dx_dt = a * x + y * z
dy_dt = b * (y-z)
dz_dt = -x*y+c*y-z
return [dx_dt, dy_dt, dz_dt]
常数是"a","b"和"c".我目前有一个每个时间步的"a"列表,我想在每个时间步骤插入,当使用scipy ode求解器时...这可能吗？
谢谢!
在不使用质量矩阵的情况下，ode 求解器ode45可以求解 y'=f(t,y)。
但是对于涉及“质量”矩阵的问题，在 ode 求解器中有一个质量矩阵选项，M(t,y)y'=f(t,y)。
“质量”矩阵究竟是什么？这个术语是否来自质量弹簧系统振荡的质量？我在文档中找不到关于此的示例代码。此外，似乎我可以在 y'=f(t,y) 等式中的 f(t,y) 中对有关 t 和 y 的信息进行编码。在什么情况/示例中会出现 M(t,y)y'=f(t,y) 需要 M(t,y) 的地方？
我正在尝试使用scipy.integrate.ode或scipy.integrate.odeint求解大量（超过一千个）初始条件的 ODE 系统，但是执行循环的速度非常慢，并且 scipy 似乎没有提供用于输入 2D 数组（由一组指定初始条件的一维数组），并且选项vectorized似乎scipy.integrate.solve_ivp并不意味着它接受初始条件的二维数组（https://docs.scipy.org/doc/scipy/reference/ generated/scipy.integrate .solve_ivp.html）。我读过一个线程询问类似的问题（Vectorized SciPy odesolver），其中一个答案建议使用scipy.integrate.odeint，但是它似乎也不接受多维数组，所以它根本不明白如何实现它。有什么解决方案可以加快进程吗？除了矢量化之外，我还考虑过使用并行计算技术，但我对此并不熟悉，而且我认为它并没有真正像矢量化那样显着地加速程序？