标签: ode

与Mathematica 8.0 NDSolve相比,使用GSL库用C语言编写的ODE求解器是否具有明显的速度优势？在准确性方面如何公平？

我的理解是编译后的代码原则上可以更快,但是现在NDSolve已经以某种方式使用了很多已编译的代码本身？

还有使用MathLink或Mathematica的编译功能来加速解决ODE的选项吗？

我正在 Matlab 中试验 ode45。我已经学会了如何将参数传递给 ode 函数，但我仍然有一个问题。假设我想计算一辆汽车的轨迹（速度曲线），并且我有一个函数，例如getAcceleration，它可以为我提供汽车的加速度以及正确的档位：[acceleration, gear] = getAcceleration(speed,modelStructure)其中modelStructure代表汽车的型号。

ode 函数将是：

function [dy] = car(t,y,modelStructure)

dy           = zeros(2,1);
dy(1)        = y(2);
[dy(2),gear] = getAcceleration(y(1),modelStructure);

然后我这样调用Ode45积分器：

tInit = 0;
tEnd  = 5,
[t,y] = ode45(@car,[tInit tEnd], [speedInitial,accelerationInitial],options,modelStructure);

问题是：如何获得矢量存储齿轮？我应该有类似的东西[t,y,gear]=ode45(....)还是应该gear在y向量内？

我一直在处理我的代码并使用事件功能，我现在能够获得汽车“齿轮”的变化（作为事件）。现在我有一个与相同代码相关的新问题。想象一下，当我评估 de 'dy' 向量时，我能够得到一个进一步的值 Z，这让我可以大幅加快调用加速度计算 (getAcceleration)：

function [dy] = car(t,y,modelStructure)

dy           = zeros(2,1);
dy(1)        = y(2);
[dy(2),Z(t)] = getAcceleration(y(1),modelStructure,Z(t-1)); 

并假设我也能够在初始条件下计算 Z。问题是我无法计算 Z 导数。

有没有办法通过 Z 值抛出步进而不集成它？

谢谢你们。

当 ODE45 的解发散时（无论原因和方式），将显示以下警告，并且求解器无法继续：

警告：t=8.190397e+01 时失败。如果不将步长减小到时间 t 时允许的最小值 (2.273737e-13) 以下，则无法满足积分容差。

我在矩阵（大量输入）上运行 ode45，所以我想自动找出哪些输入发生上述条件（失败）。我的意思是，是否有 ode45 返回的这种条件的任何其他标志可以自动写入数组？可以在if语句中使用的内容如下：

if {some variable is returned/is equal to ...} then {the solver has failed}

自动识别那些错误的输入，而无需寻找显示的警告。

我想使用deSolve R包中的显式Runge-Kutta方法ode45(别名rk45dp7)来解决具有可变步长的ODE问题.

根据deSolve文档,可以使用ode45方法而不是等距时间步长为rk求解器函数使用自适应或可变时间步长,但我不知道如何做到这一点.

rk函数被调用如下:

rk(y, times, func, parms, rtol = 1e-6, atol = 1e-6, verbose = FALSE, tcrit = NULL,
hmin = 0, hmax = NULL, hini = hmax, ynames = TRUE, method = rkMethod("rk45dp7", ... ), 
maxsteps = 5000, dllname = NULL, initfunc = dllname, initpar = parms, rpar = NULL, 
ipar = NULL, nout = 0, outnames = NULL, forcings = NULL, initforc = NULL, fcontrol = 
NULL, events = NULL, ...)

与 …

我正试图解决和涉及向量的颂歌,并且不能提出可行的答案.所以我将它分成6个分量,一个用于组件的每个时间导数,一个用于速度分量的每个时间导数.第一个值似乎是合理的,然后它跳到数百万的数字,我不知道为什么.老实说,我真的不确定如何做到这一点,我现在只是试一试.我似乎无法在网上找到任何信息,如果有这类问题的例子,可以使用一些帮助或一些链接.任何信息将非常感谢如何解决这个问题.

def dr_dt(y, t):
"""Integration of the governing vector differential equation.
d2r_dt2 = -(mu/R^3)*r with d2r_dt2 and r as vecotrs.
Initial position and velocity are given.
y[0:2] = position components
y[3:] = velocity components"""

    G = 6.672*(10**-11)
    M = 5.972*(10**24)
    mu = G*M
    r = np.sqrt(y[0]**2 + y[1]**2 + y[2]**2) 

    dy0 = y[3]
    dy1 = y[4]
    dy2 = y[5]
    dy3 = -(mu / (r**3)) * y[0]
    dy4 = -(mu / (r**3)) * y[1]
    dy5 = -(mu / (r**3)) * y[2] …

我有一个朱莉娅代码：

using DifferentialEquations\nusing Plots\nusing ParameterizedFunctions\nplotly()\nlorenz = @ode_def Lorenz begin\n  dx = \xcf\x83*(y-x)\n  dy = \xcf\x81*x-y-x*z\n  dz = x*y-\xce\xb2*z\nend \xcf\x83 = 10. \xce\xb2 = 8./3. \xcf\x81 => 28.\nu0 = [1., 5., 10.]\ntspan = (0., 2.)\nprob = ODEProblem(lorenz, u0, tspan)\nsol = solve(prob,save_timeseries=true)\nplot(sol,vars=(:x,:y,:z))\n

其结果是：\n下一个图\n
\n我如何为该图制作动画，以便它可以从 REPL 和 jupyter 运行？

我正在使用scipy食谱中的Zombie Apocalypse 示例，以了解有关使用python解决ODE系统的信息。

在此模型中，有一个方程式可以根据出生率，死亡率和初始人口提供每天的人口。然后根据人口数量，计算出有多少僵尸被制造和杀死。

我有兴趣用一系列可以告诉我们每个时间步长人口的数据列表来代替人口微分方程。我收到以下错误：

TypeError: can't multiply sequence by non-int of type 'float'

正如人们指出的那样，这是因为将单个数字乘以列表没有意义。我不确定如何在每次T时从列表中为微分方程提供数字。

这是两次尝试的代码

# solve the system dy/dt = f(y, t)
def f(y, t):
    Si = [345, 299, 933, 444, 265, 322] # replaced an equation with list
    Zi = y[0]
    Ri = y[1]
    # the model equations (see Munz et al. 2009)
    f0 = B*Si*Zi + G*Ri - A*Si*Zi
    f1 = d*Si + A*Si*Zi - G*Ri
    return [f0, f1]

我也尝试过

numbers = [345, 299, 933, …

我正在尝试对两个非线性 ODE 的系统进行数值求解。我正在使用 Scipy 的odeint功能。odeint需要一个参数y0来指定初始条件。然而，它似乎假设初始条件y0在同一时间点开始（即两个条件都在 t=0）。就我而言，我想指定为不同时间指定的两个不同边界条件（即 omega(t=0) = 0, theta(t=100) = 0）。我似乎无法弄清楚如何做到这一点，非常感谢任何帮助！

下面的一些示例代码：

from scipy.integrate import odeint

def pend(y, t, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt

b = 0.25
c = 5.0
t = np.linspace(0, 100, 101)

# I want to make these initial conditions specified at different times
y0 = [0, 0]

sol = odeint(pend, y0, t, args=(b, c))

尝试实现自适应步长 Runge-Kutta Cash-Karp 但失败并出现此错误：

home/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:15: RuntimeWarning: divide by zero encountered in double_scalars from ipykernel import kernelapp as app

我试图解决的 ODE（并在下面的示例中使用，从高阶转换为一阶 ODE 系统）如下：

这是我正在使用的实现：

def rkck(f, x, y, h, tol):
    #xn = x + h
    err = 2 * tol
    while (err > tol):
        xn = x + h
        k1 = h*f(x,y)
        k2 = h*f(x+(1/5)*h,y+((1/5)*k1)) 
        k3 = h*f(x+(3/10)*h,y+((3/40)*k1)+((9/40)*k2))
        k4 = h*f(x+(3/5)*h,y+((3/10)*k1)-((9/10)*k2)+((6/5)*k3))
        k5 = h*f(x+(1/1)*h,y-((11/54)*k1)+((5/2)*k2)-((70/27)*k3)+((35/27)*k4))
        k6 = h*f(x+(7/8)*h,y+((1631/55296)*k1)+((175/512)*k2)+((575/13824)*k3)+((44275/110592)*k4)+((253/4096)*k5))
        yn4 = y + ((37/378)*k1)+((250/621)*k3)+((125/594)*k4)+((512/1771)*k6)
        yn5 = y + ((2825/27648)*k1)+((18575/48384)*k3)+((13525/55296)*k4)+((277/14336)*k5)+((1/4)*k6)
        err = yn4[-1]-yn5[-1]
        if …

我正在尝试使用 Julia 中的 DifferentialEquations 解决谐振子。IE：

using DifferentialEquations
using Plots

m = 1.0                          
? = 1.0                     

function mass_system!(ddu,du,u,p,t)
    # a(t) = (1/m) w^2 x 
    ddu[1] = (1/m)*(?^2)*u[1]
end

v0 = 0.0                     
u0 = 1.0                  
tspan = (0.0,10.0)               

prob = SecondOrderODEProblem{isinplace}(mass_system!,v0,u0,tspan,callback=CallbackSet())
sol = solve(prob)

但它似乎并不理解 ODE 构造函数。运行后，我得到：

ERROR: LoadError: TypeError: non-boolean (typeof(isinplace)) used in boolean context
Stacktrace:
 [1] #_#219(::Base.Iterators.Pairs{Symbol,CallbackSet{Tuple{},Tuple{}},Tuple{Symbol},NamedTuple{(:callback,),Tuple{CallbackSet{Tuple{},Tuple{}}}}}, ::Type{SecondOrderODEProblem{DiffEqBas
e.isinplace}}, ::Function, ::Float64, ::Float64, ::Tuple{Float64,Float64}, ::DiffEqBase.NullParameters) at /Users/brandonmanley/.julia/packages/DiffEqBase/avuk1/src/problems/ode_problems
.jl:144
 [2] Type at ./none:0 [inlined] (repeats 2 times)
 [3] top-level scope at /Users/brandonmanley/Desktop/nBody/nBodyNN/test.jl:25 …