Building and solving a cascade of ODEs using DifferentialEquations

Specific Domains Modelling & Simulations
1

I want to build and solve a cascade of ODEs, by which term I mean that the first derivative of the k-th component of the vector x of (state) variables with respect to time t depends on x_{k}(t) itself and on the preceding x_{k-1}(t):

d/dt x_{1}(t) = f(x_{0},x_{1}) 	
d/dt x_{2}(t) = f(      x_{1},x_{2})
d/dt x_{3}(t) = f(            x_{2},x_{3})
.
.
.
d/dt x_{n}(t) = f(                           x_{n-1},x_{n}).

where x_{0}(t) is prescribed for all t in [0,tmax], that is, it is not solved for. And indeed, the f() function is identical for all the equations.

Although I am able to create such a model manually for a small number n of equations following the syntax of DifferentialEquations (see the MWE below), I was wondering if there is a more elegant (Julian) way of assembling such set of differential equations. Preferably one that also helps consequent numerical solution - perhaps the structure could be exploited somehow for solving the ODEs.

using DifferentialEquations
using Plots
pyplot()

f(x,u) = -x + u         # but could be nonlinear in general

function cascade(dx,x,x0,t)
 dx[1] = f(x[1],x0(t))
 dx[2] = f(x[2],x[1])
 dx[3] = f(x[3],x[2])
 dx[4] = f(x[4],x[3])
 dx[5] = f(x[5],x[4])
end

xinit = zeros(5)
ti = 0.0
tf = 10.0
tspan = (ti,tf)

x₀(t) = 0.5-t >= 0.0 ? 0.0 : 1.0

prob = ODEProblem(cascade,xinit,tspan,x₀)
sol = solve(prob)

plot(sol.t,x₀,linewidth=2,label="x0")
plot!(sol,linewidth=2,xaxis="t",yaxis="x(t)",label=["x1" "x2" "x3" "x4" "x5"])
2

This works fine (not sure how to exploit structure though):

function cascade(dx, x, x0, t)
    dx[1] = f(x[1], x0(t))
    for i in 2:length(x)
        dx[i] = f(x[i], x[i-1])
    end
end
2 Likes
3

Thanks.

What motivates me to think about exploiting the structure of the problem is that in order to solve for x_i(), only x_{i-1}() needs to be known (as an input or parameter for that matter), hence a first-order ODE. The whole problem could be thus solved as a sequence of first-order nonhomogeneous ordinary differential equations. The current approach builds and solves a (possibly) huge monotitic system of ODEs.

I am certainly able to hard-wire it somehow, I was just wondering if somebody can immediately recognize this as an opportunity for some advanced Julia technique (perhaps designing some iterators?). Just a hint, I will then explore it on my own.

4

The Jacobian has a special structure that you can pass to DifferentialEquations

1 Like
5

If you’re looking to take advantage of the structure for the sake of terseness, the accumulate! function will get you there. It sacrifices some readability though, in my opinion, because you have to do a little extra mental accounting.

If you’re looking to take advantage of the structure for the sake of performance, as @rveltz said, the Jacobian can be defined explicitly and passed into ODEFunction. Here is a quick example of both using x₀(t) = 3*sqrt(t) as the boundary function and f(x1, x2) = -x1^2 * sin(x2) as the cascading function:

using DifferentialEquations

tspan = (0.0, 1.0)
xinit = ones(5)
x₀(t) = 3*sqrt(t)

f(x1, x2) =        -x1^2 * sin(x2)
df_dx1(x1, x2) = -2*x1   * sin(x2)
df_dx2(x1, x2) =   -x1^2 * cos(x2)

cascade!(dx, x, x0, t) = accumulate!(f, dx, x; init=x0(t))

function jac(J, x, x0, t)
    J[1,1] = df_dx1(x[1], x0(t))
    for i in 2:length(x)
        J[i,   i] = df_dx1(x[i], x[i-1])
        J[i, i-1] = df_dx2(x[i], x[i-1])
    end
end

ff = ODEFunction(cascade!; jac=jac)

prob = ODEProblem(ff, xinit, tspan, x₀)
2 Likes
6

And then you can also pass the jac_prototype for the sparsity pattern, which will be bidiagonal.

3 Likes