我目前正在尝试使用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有任何建议？
我是Matlab的新手.我希望你能帮助我.我必须使用ODE45函数来解决一个ODE系统.这是描述我的等同的功能.
function dNdt = rateEquations(t, y)
%populations of corresponding state
Ng = y(1);
Ns = y(2);
Nt = y(3);
%All constants used are dropped for the sake of easy reading.
注意参数F.
%rate equations
dNs = s0 * Ng * F - Ns/ t_S1;
dNt = Ns / t_ISC - Nt / t_T1;
dNg = -dNt - dNs;
dNdt = [dNg; dNs; dNt];
end
然后,在我的脚本.m文件中,我在'for循环'中调用ode45函数.在每次迭代期间,我必须更改参数F并将其传递给我的'rateEquations' - 函数.但我不知道如何实现它.
for T = Tmin: dt : Tmax
%initial conditions
initialConditions = [N0 0 0]; …我有以下功能包含一些颂歌:
myfunction <- function(t, state, parameters) {
with(as.list(c(state, parameters)),{
if (X>20) { # this is an internal threshold!
Y <- 35000
dY <- 0
}else{
dY <- b * (Y-Z)
}
dX <- a*X^6 + Y*Z
dZ <- -X*Y + c*Y - Z
# return the rate of change
list(c(dX, dY, dZ),Y,dY)
})
}
以下是一些结果:
library(deSolve)
parameters <- c(a = -8/3, b = -10, c = 28)
state <- c(X = 1, Y = 1, Z = 1)
times <- …我有一个ordinary differential equation及时的域名如下:
C*du/dt = -g*u + I
哪里 I = A*t/tau*exp^(1-t/tau)
在freq域中:
u(w) = I(w)/(g*(1+C/g*j*w))
j 是一个复杂的数字 sqrt(-1)
因此,我可以u(t)通过快速傅里叶变换(fft)进入频域 ,然后再使用ifft.
代码:
t = np.linspace(0.,499.9,5000)
I = q*t*np.exp(1-t/Tau_ca)/Tau_ca
u1 = np.fft.ifft(np.fft.fft(I)/(G_l*(1.+1.j*(C_m/G_l)*np.fft.fftfreq(t.shape[-1]))))
然而,当我将其u(t)与其他方法(如微分方程的数值积分或其分析形式)进行比较时,这是不正确的.我已经尝试过并且未能成功找出我的错误所在.
请指教.
y''' + 41y'' + 360y' + 900y = 600x' + 1200x;
y(0)= 2 ; y'(0)= 1 ; y''(0) = -0.05
如何使用ODE45函数求解该等式？
我试过这个:
==>
function dydt=f(t,y)
dydt = [y(2) ; y(3) ; -41*y(3)-360*y(2)- 900*y(1)]
==>
clear all;
timerange=[0 1.4]; %seconds
initialvalues=[2 1 -0.05];
[t,y]=ode45(@dydt, timerange, initialvalues)
plot(t,y(:,1));
但我需要将X部分放在等式中 - 我不知道如何......
使用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 …我有一个简单的二阶ODE的同类解决方案,当我尝试使用Sympy求解初始值时,返回相同的解决方案.它应该替代y(0)和y'(0)并产生没有常数的解,但不能.这是设置方程的代码(它是弹簧平衡方程,k =弹簧常数,m =质量).我在其他地方使用的一些冗余符号,抱歉.
%matplotlib inline
from sympy import *
m,k, x,t,s,T, omega,A,B = symbols('m k x t s T omega A B',float=True)
a = symbols('a', positive=True)
f,y,g,H,delta,S=symbols('f y g H delta S',cls=Function)
Eq1 = Eq(m*diff(y(t),t,2)+k*y(t))
Eq1
结果是(正确):$ y {\ left(t\right)} = C_ {1} e ^ { - t\sqrt { - \frac {k} {m}}} + C_ {2} e ^ { t\sqrt { - \frac {k} {m}}} $
y(t)= C1e ^( - t√( - k/m))+ C2e ^(t√( - km)),也有y_n = c1.cos(√( - …
我想从数字的角度理解混合建模(特别是状态事件)背后的一般思路(虽然我不是数学家:)).鉴于以下Modelica模型:
model BouncingBall
constant Real g=9.81
Real h(start=1);
Real v(start=0);
equation
der(h)=v;
der(v)=-g;
algorithm
when h < 0 then
reinit(v,-pre(v));
end when;
end BouncingBall;
我理解的概念when和reinit.
when声明中的等式仅在条件成立时才有效吗？2sec.由于我使用多步求解器,这意味着求解器"超过2秒",识别出h<0(假设在模拟time = 2.5sec中h = -0.7).这是什么意思"使用交叉函数搜索事件的时间？是否有简单的解释(示例)？pre()在这种情况下,操作意味着什么？noEvent():"表达式是字面意义而不是生成交叉函数.由于没有交叉函数,因此没有要求表达式可以在事件限制之外进行评估":这是什么意思？给出与弹跳球相同的例子:求解器在时间2.5检测到h = 0.7.有没有区别noEvent()？我试图使用nlme和拟合一阶微分模型lsoda。这是基本思想：我首先定义允许生成微分方程解的函数：
library(deSolve)
ODE1 <- function(time, x, parms) {with(as.list(c(parms, x)), {
import <- excfunc(time)
dS <- import*k/tau - (S-yo)/tau
res <- c(dS)
list(res)})}
solution_ODE1 = function(tau1,k1,yo1,excitation,time){
excfunc <- approxfun(time, excitation, rule = 2)
parms <- c(tau = tau1, k = k1, yo = yo1, excfunc = excfunc)
xstart = c(S = yo1)
out <- lsoda(xstart, time, ODE1, parms)
return(out[,2])
}
然后，根据两个ID的公式生成数据：
time <- 0:49
excitation <- c(rep(0,10),rep(1,10),rep(0,10),rep(1,10),rep(0,10))
simu_data <- data.frame(signal = c(solution_ODE1(3,2,0.1,excitation,time)+rnorm(length(time),0,0.1),
solution_ODE1(3.2,1.5,0.3,excitation,time)+rnorm(length(time),0,0.1)),
time = rep(time,2),
excitation …我正在尝试使用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，但是它似乎也不接受多维数组，所以它根本不明白如何实现它。有什么解决方案可以加快进程吗？除了矢量化之外，我还考虑过使用并行计算技术，但我对此并不熟悉，而且我认为它并没有真正像矢量化那样显着地加速程序？