Решить систему ОДУ с чтением во внешнем форсировании

Question

В Юлии я хочу решить систему ОДУ с внешним принуждением g1(t), g2(t) лайк

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

с принуждениями, считанными из файла.

Я использую это исследование для изучения Юлии и пакета diffrentialEquations, но мне трудно найти правильный подход.

Я мог бы представить, что с помощью 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))