# How to solve a combination of ODEs and PDEs and organize the solution

**URL:** <https://discourse.julialang.org/t/how-to-solve-a-combination-of-odes-and-pdes-and-organize-the-solution/95092>\
**Category:** General Usage\
**Tags:** question, differentialequation\
**Created:** [February 23, 2023, 4:55pm UTC](https://discourse.julialang.org/t/how-to-solve-a-combination-of-odes-and-pdes-and-organize-the-solution/95092 "2023-02-23T16:55:40Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![markowkes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/markowkes/32/43198_2.png) [@markowkes](https://discourse.julialang.org/u/markowkes)\
**Post date:** [February 23, 2023, 4:55pm UTC](https://discourse.julialang.org/t/how-to-solve-a-combination-of-odes-and-pdes-and-organize-the-solution/95092/1 "2023-02-23T16:55:40Z")

</div>

Is it possible to use a struct to hold the dependent variables when using the solvers in DifferentialEquations.jl, or is it needed to create a single vector u of all dependent variables?

For example, I have a system with 3 ODEs

- dXt/dt = …
- dSt/dt = …
- dLf/dt = …

and 2 PDES that depend on one spatial dimension

- dPb(z)/dt = …
- dSb(z)/dt = …

I’m solving these by discretizing the spatial derivatives in the PDEs by creating a grid with N grid cells and using finite differences to create N ODEs representing each PDE. This gives a total of 3 + 2\*N ODEs that are organized into a single vector solved with DifferentialEquations.jl. Everything works, but the solution is difficult to work with since it requires splitting the variables into Xt, St, Lf, Pb, and Sb. I should note that it’s actually a little more complicated, but this is the general idea.

It would be much easier to have a struct that organizes the solution, e.g.,

```julia
struct sol 
    Xt :: Float
    St :: Float
    Lf :: Float
    Pb :: Vector{Float}
    Sb :: Vector{Float}
end

```

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [February 23, 2023, 5:01pm UTC](https://discourse.julialang.org/t/how-to-solve-a-combination-of-odes-and-pdes-and-organize-the-solution/95092/2 "2023-02-23T17:01:00Z")

</div>

I have used [ComponentArrays](https://github.com/jonniedie/ComponentArrays.jl) to do similar things, I only used an ODE solver but I suspect it works with the PDE solvers as well.

---

<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:** [February 23, 2023, 5:02pm UTC](https://discourse.julialang.org/t/how-to-solve-a-combination-of-odes-and-pdes-and-organize-the-solution/95092/3 "2023-02-23T17:02:46Z")

</div>

> [@markowkes](#):
>
> Is it possible to use a struct to hold the dependent variables when using the solvers in DifferentialEquations.jl, or is it needed to create a single vector u of all dependent variables?

I usually just use a bigger vector, like is shown in the showcase example:

[https://docs.sciml.ai/Overview/stable/showcase/gpu\_spde/](https://docs.sciml.ai/Overview/stable/showcase/gpu_spde/)

However, a great library for writing this kind of code is ComponentArrays.jl which then allows the pieces to be named.

> **[GitHub - jonniedie/ComponentArrays.jl: Arrays with arbitrarily nested named...](https://github.com/jonniedie/ComponentArrays.jl)**
>
> Arrays with arbitrarily nested named components. Contribute to jonniedie/ComponentArrays.jl development by creating an account on GitHub.

Writing your PDE+ODE with component arrays is pretty natural and gives an efficient form.

---

<div class="post-metadata">

**Author:** ![markowkes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/markowkes/32/43198_2.png) [@markowkes](https://discourse.julialang.org/u/markowkes)\
**Post date:** [February 24, 2023, 3:47pm UTC](https://discourse.julialang.org/t/how-to-solve-a-combination-of-odes-and-pdes-and-organize-the-solution/95092/4 "2023-02-24T15:47:27Z")

</div>

Thanks @contradict and @ChrisRackauckas for the suggestion to use [ComponentArrays](https://github.com/jonniedie/ComponentArrays.jl). I have implemented it in my code and it works well.

I did find it non-trival to access the parts of the DiffEq solution using the component labels. Ultimately, I found the `label2index()` function which is part of ComponentArrays can be used to find the appropriate indices which can be used with the provided solution analysis tools or plot recipes in DifferentialEquations.jl.

If there is a way to access the components without using `label2index()`, please share!

Below is a MWE that highlights how this can work.

```julia

using ComponentArrays
using OrdinaryDiffEq
using Parameters: @unpack
using Plots

# Define function to get index of 
# each labeled component in solution
idx(sol,var) = label2index(sol.u[1],var)

function odetest()
    # Define and solve 
    # ODE dXo/dt = Xo and 
    # PDE d(Xp=A,B,C)/dt = Xp
    # ---------------------

    # Define right-hand-side of ODEs
    function dsol_dt!(dsol, sol, p, t)
        @unpack Xo,Xp = sol
        dsol.Xo = Xo
        dsol.Xp = Xp
        return nothing
    end

    # Initial conditions
    Xo₀ = 1.0 # 
    Xp₀ = ComponentArray(A = 1.0, B = 2.0, C = 3.0)
    sol₀ = ComponentArray(Xo=Xo₀,Xp=Xp₀)

    # Setup problem and solve
    tspan = (0.,1.0)
    prob = ODEProblem(dsol_dt!,sol₀,tspan)
    sol = solve(prob,Tsit5(),)

    # Access final solution 
    # --------------------- 
    
    # Access components
    Xo = sol[idx(sol,"Xo"),:]
    Xp = sol[idx(sol,"Xp"),:]

    # Access components with interpolation
    t = 0.0:0.1:1.0
    Xo_dense = sol(t,idxs = idx(sol,"Xo"))
    Xp_dense = sol(t,idxs = idx(sol,"Xp"))

    # Plot solution and components
    p1 = plot(sol ,title="Full solution")
    p2 = plot(sol,idxs = idx(sol,"Xo"),title="Only Xo")
    p3 = plot(sol,idxs = idx(sol,"Xp"),title="Only Xp")
    p = plot(p1,p2,p3,layout=(1,3))
    display(p)
end
odetest()

```

![odetest](https://global.discourse-cdn.com/julialang/original/3X/6/4/64ca6c1332be7c449b4ed008c3f775865fec10f4.png)

---

<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:** [February 25, 2023, 3:23pm UTC](https://discourse.julialang.org/t/how-to-solve-a-combination-of-odes-and-pdes-and-organize-the-solution/95092/5 "2023-02-25T15:23:20Z")

</div>

> [@markowkes](#):
>
> `Xo = sol[idx(sol,"Xo"),:]`

`sol[end].Xo`?

---

<div class="post-metadata">

**Author:** ![markowkes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/markowkes/32/43198_2.png) [@markowkes](https://discourse.julialang.org/u/markowkes)\
**Post date:** [February 26, 2023, 12:22am UTC](https://discourse.julialang.org/t/how-to-solve-a-combination-of-odes-and-pdes-and-organize-the-solution/95092/6 "2023-02-26T00:22:47Z")

</div>

That does work to access the solution at the final time, but access more than one time, e.g. at all the times with `Xp = sol[1:end].Xp`, does not work. This gives the error `ERROR: type ODESolution has no field Xp`

---

<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:** [February 26, 2023, 1:38am UTC](https://discourse.julialang.org/t/how-to-solve-a-combination-of-odes-and-pdes-and-organize-the-solution/95092/7 "2023-02-26T01:38:57Z")

</div>

It works on each i, so just use a compeehension