# Building and solving a cascade of ODEs using DifferentialEquations

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

dx = f(x,x0(t))
dx = f(x,x)
dx = f(x,x)
dx = f(x,x)
dx = f(x,x)
end

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

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

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"])
``````

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

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

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.

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

1 Like

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

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