通过读入外力求解 ODE 系统

Question

在 Julia 中，我想求解具有外力的 ODE 系统 g1(t), g2(t)，例如

dx1(t) / dt = f1(x1, t) + g1(t)
dx2(t) / dt = f2(x1, x2, t) + g2(t)

使用从文件中读入的强迫。

我正在使用这项研究来学习 Julia 和微分方程包，但我很难找到正确的方法。

我可以想象使用 callback 可以工作，但这似乎很麻烦。

你知道如何实现这样的外力吗？

Answer 1

您可以在集成函数中使用函数。因此，您可以使用 Interpolations.jl 之类的方法从文件中的数据构建插值多项式，然后执行以下操作：

g1 = interpolate(data1, options...)
g2 = interpolate(data2, options...)
p = (g1,g2) # Localize these as parameters to the model

function f(du,u,p,t)
g1,g2 = p
  du[1] = ... + g1[t] # Interpolations.jl interpolates via []
  du[2] = ... + g2[t]
end
# Define u0 and tspan
ODEProblem(f,u0,tspan,p)

Answer 2

感谢@Chris Rackauckas 提出的一个很好的问题和很好的回答。下面是此类问题的完整可重现示例。请注意，Interpolations.jl 已将索引更改为 g1(t)。

using Interpolations
using DifferentialEquations
using Plots

time_forcing = -1.:9.
data_forcing = [1,0,0,1,1,0,2,0,1, 0, 1]
g1_cst = interpolate((time_forcing, ), data_forcing, Gridded(Constant()))
g1_lin = scale(interpolate(data_forcing, BSpline(Linear())), time_forcing)

p_cst = (g1_cst) # Localize these as parameters to the model
p_lin = (g1_lin) # Localize these as parameters to the model


function f(du,u,p,t)
  g1 = p
  du[1] = -0.5 + g1(t) # Interpolations.jl interpolates via ()
end

# Define u0 and tspan
u0 = [0.]
tspan = (-1.,9.) # Note, that we would need to extrapolate beyond 
ode_cst = ODEProblem(f,u0,tspan,p_cst)
ode_lin = ODEProblem(f,u0,tspan,p_lin)

# Solve and plot
sol_cst = solve(ode_cst)
sol_lin = solve(ode_lin)

# Plot
time_dense = -1.:0.1:9.
scatter(time_forcing, data_forcing,   label = "discrete forcing")
plot!(time_dense, g1_cst(time_dense), label = "forcing1",  line = (:dot,   :red))
plot!(sol_cst,                        label = "solution1", line = (:solid, :red))
plot!(time_dense, g1_lin(time_dense), label = "forcing2",  line = (:dot,   :blue))
plot!(sol_lin,                        label = "solution2", line = (:solid, :blue))

通过读入外力求解 ODE 系统

solve system of ODEs with read in external forcing

julia

ode

differentialequations.jl