solve system of ODEs with read in external forcing

Question

In Julia, I want to solve a system of ODEs with external forcings g1(t), g2(t) like

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

with the forcings read in from a file.

I am using this study to learn Julia and the package DifferentialEquations, but I am having difficulties finding the correct approach.

I could imagine that using a callback could work, but that seems pretty cumbersome.

Do you have an idea of how to implement such an external forcing?

score 5 · Accepted Answer · answered Mar 22 '18 at 13:04

5

You can use functions inside of the integration function. So you can use something like Interpolations.jl to build an interpolating polynomial from the data in your file, and then do something like:

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)

answered Mar 22 '18 at 13:04

Chris Rackauckas

18,645
3
50
81

It's hard to know without seeing what you're trying to do. My guess would be you may need to re-write the data into a different format: maybe it only wants to interpolate vectors. We can discuss this in more detail in the [Gitter channel](https://gitter.im/JuliaLang/julia) the [Julia Slack](https://slackinvite.julialang.org/), or the [Julia Discourse](https://discourse.julialang.org/). These avenues are better for long discussions, whereas StackOverflow is just simple question->answer. – Chris Rackauckas Mar 22 '18 at 16:56
OK, thanks, I had some problems with reloading the new code. I'll have to play around a bit more and in case I'll stumble upon new problems, I'll had over to the discussion groups. Thanks again! – LSchueler Mar 22 '18 at 17:01

score 2 · Answer 2 · answered May 07 '20 at 12:50

Thanks for a nice question and nice answer by @Chris Rackauckas. Below a complete reproducible example of such a problem. Note that Interpolations.jl has changed the indexing to 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))

solve system of ODEs with read in external forcing

2 Answers2

Linked