1

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?

2 Answers 2

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)
Sign up to request clarification or add additional context in comments.

2 Comments

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 the Julia Slack, or the Julia Discourse. These avenues are better for long discussions, whereas StackOverflow is just simple question->answer.
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!
2

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))

enter image description here

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.