# Question on ODE system with inhomogeneous coefficients, parameters, code design

**URL:** <https://discourse.julialang.org/t/question-on-ode-system-with-inhomogeneous-coefficients-parameters-code-design/56849>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, forwarddiff\
**Created:** [March 10, 2021, 4:04am UTC](https://discourse.julialang.org/t/question-on-ode-system-with-inhomogeneous-coefficients-parameters-code-design/56849 "2021-03-10T04:04:43Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ininterrompue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ininterrompue/32/5594_2.png) [@Ininterrompue](https://discourse.julialang.org/u/Ininterrompue)\
**Post date:** [March 10, 2021, 4:04am UTC](https://discourse.julialang.org/t/question-on-ode-system-with-inhomogeneous-coefficients-parameters-code-design/56849/1 "2021-03-10T04:04:43Z")

</div>

Hi,

I have a question relating to code design. I’m using some ODE solvers to calculate a system of ODEs with time-varying coefficients which I am interested in. They take on the form `F`, `G`, and `dG` in the MWE and depend on some parameters which I’ve defined in a separate struct. The fact that the functions are listed twice seems to be an obvious case of code smell, but there’s a couple issues

- `ForwardDiff.derivative` only takes unary arguments, but e.g. `G` depends on parameters `b` and `c` (the actual function also depends on something like `u[1] + u[2]`)
- The actual function to differentiate depends on `t` in a much more complicated manner than `G` below, with several layers of functions in `t`. `G` is needed in the solver, but once the solver is done, I don’t have access to `G` without copying the function definition.

In the real code, there are about a dozen parameters and a dozen functions like the three I’ve given, so I’d like to only write them out once. Is this a case of using the integrator interface? Any assistance is greatly appreciated.

```julia
using OrdinaryDiffEq, Plots
import ForwardDiff

struct Parameters{T}
    a::T
    b::T
    c::T
end

function f(du, u, p, t)
    F(t) = sqrt(p.a * t)
    G(t) = p.b - exp(-p.c * t)
    dG(t) = ForwardDiff.derivative(G, t)

    du[1] = -F(t) * u[2] + dG(t)
    du[2] = +G(t) * u[1]
end

p = Parameters(0.5, 0.7, 0.2)
u0 = [1.0, 0.0]; tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob, Tsit5())

F(t) = sqrt(p.a * t)
G(t) = p.b - exp(-p.c * t)
dG(t) = ForwardDiff.derivative(G, t)

# plot F, G, dG, etc...
plt = plot(F, tspan...)
plot!(G, tspan...)
display(plt)

```

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [March 10, 2021, 4:41am UTC](https://discourse.julialang.org/t/question-on-ode-system-with-inhomogeneous-coefficients-parameters-code-design/56849/2 "2021-03-10T04:41:38Z")

</div>

An easy alternative: modify your functions to take your parameter struct as an argument:

```julia
F(p, t) = sqrt(p.a * t)
G(p, t) = p.b - exp(-p.c * t)
dG(p, t) = ForwardDiff.derivative(t -> G(p, t), t)

function f(du, u, p, t)
    du[1] = -F(p, t) * u[2] + dG(p, t)
    du[2] = G(p, t) * u[1]
end

```

---

<div class="post-metadata">

**Author:** ![Ininterrompue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ininterrompue/32/5594_2.png) [@Ininterrompue](https://discourse.julialang.org/u/Ininterrompue)\
**Post date:** [March 10, 2021, 4:48am UTC](https://discourse.julialang.org/t/question-on-ode-system-with-inhomogeneous-coefficients-parameters-code-design/56849/3 "2021-03-10T04:48:42Z")

</div>

> `dG(p, t) = ForwardDiff.derivative(t -> G(p, t), t)`

Thanks, I will keep this type of line in mind.

However I must ask what to do in the case that `G` or `dG` depends on `u` as well. For example if I multiplied `G` by `u[1] + u[2]`.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [March 10, 2021, 4:55am UTC](https://discourse.julialang.org/t/question-on-ode-system-with-inhomogeneous-coefficients-parameters-code-design/56849/4 "2021-03-10T04:55:38Z")

</div>

Just pass another argument:

```julia
G(u, p, t) = (u[1] + u[2])*(p.b - exp(-p.c * t))
dG(u, p, t) = ForwardDiff.derivative(t -> G(u, p, t), t)

```

A separate note: for small state vectors, you can get a decent speed-up by using a `StaticVector` instead of a `Vector`.

```julia
function g(u, p, t)
    du1 = -F(p, t) * u[2] + dG(p, t)
    du2 = G(p, t) * u[1]
    SA[du1, du2]
end

prob_f = ODEProblem(f, u0, tspan, (a = 0.5, b = 0.7, c = 0.2))
prob_g = ODEProblem(g, SA[1.0, 0.0], tspan, (a = 0.5, b = 0.7, c = 0.2))

```

```julia
julia> @btime solve($prob_f, Tsit5());
  8.700 μs (94 allocations: 12.17 KiB)

julia> @btime solve($prob_g, Tsit5());
  3.900 μs (26 allocations: 5.56 KiB)

```

---

<div class="post-metadata">

**Author:** ![Ininterrompue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ininterrompue/32/5594_2.png) [@Ininterrompue](https://discourse.julialang.org/u/Ininterrompue)\
**Post date:** [March 10, 2021, 6:01am UTC](https://discourse.julialang.org/t/question-on-ode-system-with-inhomogeneous-coefficients-parameters-code-design/56849/5 "2021-03-10T06:01:27Z")

</div>

Much appreciated. I’ll take a look at the static vectors; I suspect performance may become important in the future.

While this thread is fresh, can I ask how one would implement e.g. if `F` was given by an integral, something like F(t) = \int\_{0}^{t}u\_{1}^{2}(t')\,\mathrm{d}t', which requires information from the initial time up to the present time during each stop of the solve? It’s a straightforward problem if I coded the time integration manually and preallocated a bunch of vectors, but it’d be useful to know how to do it with ODEProblem.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [March 10, 2021, 6:08am UTC](https://discourse.julialang.org/t/question-on-ode-system-with-inhomogeneous-coefficients-parameters-code-design/56849/6 "2021-03-10T06:08:38Z")

</div>

The easiest thing is to include it in your state vector like `[u1, u2, F]` - you’re integrating the ODE anyways, so why not just take care of it there?