从 ODE 问题中提取导数 (du/dt)

Question

我正在解决多个 ODE 问题。我需要解 (u) 和解的导数 (du)。对于较小的 ODE，我可以执行以下操作

using DifferentialEquations


function SB(du,u,p,t)
    du[1]=@. u[2]
    du[2]=@. ((-0.5*u[2]^2)*(3-u[2]/(p[4]))+(1+(1-3*p[7])*u[2]/p[4])*((p[6]-p[5])/p[2]+2*p[1]/(p[2]*p[9]))*(p[9]/u[1])^(3*p[7])-2*p[1]/(p[2]*u[1])-4*p[3]*u[2]/(p[2]*u[1])-(1+u[2]/p[4])*(p[6]-p[5]+p[10]*sin(2*pi*p[8]*t))/p[2]-p[10]*u[1]*cos(2*pi*p[8]*t)*2*pi*p[8]/(p[2]*p[4]))/((1-u[2]/p[4])*u[1]+4*p[3]/(p[2]*p[4]))
end

R0=2e-6
ps=250e3
f=2e6


u0=([R0 0])
tspan=(0,100/f)

p=[0.0725, 998, 1e-3,1481, 0, 1.01e5,7/5,f, R0, ps]

prob = ODEProblem(SB,u0,tspan,p)


@time u = solve(prob,Tsit5(),alg_hints=[:stiff],saveat=0.01/f,reltol=1e-8,abstol=1e-8)
t=u.t
u2=@. ((-0.5*u[2,:]^2)*(3-u[2,:]/(p[4]))+(1+(1-3*p[7])*u[2,:]/p[4])*((p[6]-p[5])/p[2]+2*p[1]/(p[2]*p[9]))*(p[9]/u[1,:])^(3*p[7])-2*p[1]/(p[2]*u[1,:])-4*p[3]*u[2,:]/(p[2]*u[1,:])-(1+u[2,:]/p[4])*(p[6]-p[5]+p[10]*sin(2*pi*p[8]*t))/p[2]-p[10]*u[1,:]*cos(2*pi*p[8]*t)*2*pi*p[8]/(p[2]*p[4]))/((1-u[2,:]/p[4])*u[1,:]+4*p[3]/(p[2]*p[4]))

其中 u2 基本上是 SB 函数中的 du[2]。随着我的 ODE 大小的增长（> 500 个耦合 ODE 和 > 500X500 矩阵），这很快变得不切实际。有没有办法要求 DifferentialEquations.jl 包（或任何其他方式）在求解 ODE 时导出 du[i]s？我了解到 DiffEqSensitivity.jl 包能够提供 du/dps 来检查模型对 p 的敏感度。有没有类似extract du/dts的东西？

Answer 1

我会一起使用两个不同的组件。首先，当您使用非常大的 ODE 时，您将只想保存解决方案的特定部分或减少的部分。对于这一点，SavingCallback是非常有帮助的。

http://diffeq.sciml.ai/latest/features/callback_library#SavingCallback-1

例如，以下求解 ODE，并且仅保存每一步解的迹和范数：

using DiffEqCallbacks, OrdinaryDiffEq, LinearAlgebra
prob = ODEProblem((du,u,p,t) -> du .= u, rand(4,4), (0.0,1.0))
saved_values = SavedValues(Float64, Tuple{Float64,Float64})
cb = SavingCallback((u,t,integrator)->(tr(u),norm(u)), saved_values)
sol = solve(prob, Tsit5(), callback=cb)

现在您可以使用它来保存您需要的内容。第二部分是使用integrator来获得导数。您可以看到get_du!可用于提取当前（已计算的）导数：

http://diffeq.sciml.ai/latest/basics/integrator#Misc-1

此外，您可以利用积分器上的插值。integrator(t,Val{1})将给出当前解的一阶导数t。