# Building and solving a cascade of ODEs using DifferentialEquations

**URL:** <https://discourse.julialang.org/t/building-and-solving-a-cascade-of-odes-using-differentialequations/29963>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [October 15, 2019, 11:59pm UTC](https://discourse.julialang.org/t/building-and-solving-a-cascade-of-odes-using-differentialequations/29963 "2019-10-15T23:59:50Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [October 15, 2019, 11:59pm UTC](https://discourse.julialang.org/t/building-and-solving-a-cascade-of-odes-using-differentialequations/29963/1 "2019-10-15T23:59:50Z")

</div>

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

```nohighlight
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.

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

```

---

<div class="post-metadata">

**Author:** ![IlyaOrson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ilyaorson/32/1681_2.png) [@IlyaOrson](https://discourse.julialang.org/u/IlyaOrson)\
**Post date:** [October 16, 2019, 5:25am UTC](https://discourse.julialang.org/t/building-and-solving-a-cascade-of-odes-using-differentialequations/29963/2 "2019-10-16T05:25:39Z")

</div>

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [October 16, 2019, 11:43am UTC](https://discourse.julialang.org/t/building-and-solving-a-cascade-of-odes-using-differentialequations/29963/3 "2019-10-16T11:43:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [October 16, 2019, 1:10pm UTC](https://discourse.julialang.org/t/building-and-solving-a-cascade-of-odes-using-differentialequations/29963/4 "2019-10-16T13:10:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jonniedie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonniedie/32/12842_2.png) [@jonniedie](https://discourse.julialang.org/u/jonniedie)\
**Post date:** [October 16, 2019, 6:12pm UTC](https://discourse.julialang.org/t/building-and-solving-a-cascade-of-odes-using-differentialequations/29963/5 "2019-10-16T18:12:05Z")

</div>

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:

```julia
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₀)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 16, 2019, 6:49pm UTC](https://discourse.julialang.org/t/building-and-solving-a-cascade-of-odes-using-differentialequations/29963/6 "2019-10-16T18:49:49Z")

</div>

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