1

我有多个要解决的 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 的敏感性。有没有类似于提取 du/dts 的东西？

4

2 回答 2

5

我会一起使用两个不同的组件。首先，当您获得非常大的 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。

于 2019-07-13T10:55:30.390 回答
-2

@ChrisRackauckas 我确实需要定义求解器来解决的每个时间步骤。

get_du!(out,integrator) 给了我一个数组，其中所有点都具有相同的值。我在某个地方犯了错误吗？

prob = ODEProblem(SB,u0,tspan,p) Rdot=zeros(50001,2) u = init(prob,SSPRK22(),dt=1e-9,reltol=1e-8,abstol=1e-8) solve!(u) get_du!(Rdot,u) U=u.sol 基本上第二个输出（du [2]）的导数必须等于我在上一篇文章中定义的u2。`

于 2019-07-17T19:45:31.433 回答