Julia DifferentialEquations.jl 所有变量输出

Question

我有以下例子：

using DifferentialEquations

function test1(du,u,p,t)
    a,b,c = p
    
    d=a^0.1*(t+1)
    e=u[1]/a
    f=u[2]/d
    
    du[1] = a*u[1]
    du[2] = d*u[2]
    du[3] = b*u[2] - c*u[3]
end

p = (2,0.75,0.8)
u0 = [1.0;1.0;1.0]
tspan = (0.0,3.0)
prob = ODEProblem(test1,u0,tspan,p)
sol = solve(prob,saveat=0.3)

sol

对象包含状态输出，但是，我还需要有效的其他变量（“d”，“e”，“f”）。

我能得到的最接近的是：

function test2(du,u,p,t)
    
    global i
    global Out_values
    global sampletimes
    
    a,b,c = p
    
    d=a^0.1*(t+1)
    e=u[1]/a
    f=u[2]/d
    
    if t in sampletimes
        
        Out_values[1,i] = d
        Out_values[2,i] = e
        Out_values[3,i] = f
        i=i+1
    
    end
    
    du[1] = a*u[1]
    du[2] = d*u[2]
    du[3] = b*u[2] - c*u[3]
end

sampletimes = tspan[1]:0.3:tspan[2]
Out_values = Array{Float64}(undef, 3, 2*11)
i=1

prob = ODEProblem(test2,u0,tspan,p)
sol = solve(prob,saveat=0.3,tstops=sampletimes)

但是，这个解决方案并不理想，因为：

1. 它重复

   saveat

   ，我得到两组略有不同的输出（不知道为什么），并且
2. 如果我决定不使用

   saveat

   并且我想输出所有解决方案，即

   sol = solve(prob)

   ，它就无法扩展。

如有任何帮助，我们将不胜感激。

Answer 1

如果您以 ModelingToolkit 形式编写此模型：

using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D

@variables 
@parameters 

@mtkmodel SampleODE begin
  @parameters begin
    a=2 
    b=0.75 
    c=0.8
  end
  @variables begin
    d(t) 
    e(t) 
    f(t) 
    x(t)=1.0 
    y(t)=1.0 
    z(t)=1.0
  end
  @equations begin
    d ~ a^0.1*(t+1)
    e ~ x/a
    f ~ y/d
    D(x) ~ a*x
    D(y) ~ d*y
    D(z) ~ b*y - c*z
  end
end

tspan = (0.0,3.0)
@mtkbuild sys = SampleODE()
prob = ODEProblem(sys,[],tspan)
sol = solve(prob,saveat=0.3)

然后您可以使用符号接口来求解计算任何输出。例如：

julia> sol[sys.d]
11-element Vector{Float64}:
 1.0717734625362931
 1.3933055012971811
 1.714837540058069
 ⋮
 3.6440297726233966
 3.9655618113842848
 4.2870938501451725