# Solving a non-linear P(I)DE - age-structured equation

**URL:** https://discourse.julialang.org/t/solving-a-non-linear-p-i-de-age-structured-equation/132677
**Category:** Numerics
**Tags:** pde, numerics
**Created:** [September 26, 2025, 2:34pm UTC](https://discourse.julialang.org/t/solving-a-non-linear-p-i-de-age-structured-equation/132677 "2025-09-26T14:34:53Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![jucheval](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jucheval/32/213961_2.png) [@jucheval](https://discourse.julialang.org/u/jucheval)
#### Post date: [September 26, 2025, 2:34pm UTC](https://discourse.julialang.org/t/solving-a-non-linear-p-i-de-age-structured-equation/132677/1 "2025-09-26T14:34:53Z")

</div>

Hi everyone !

I am searching for some way to solve the following non-linear P(I)DE, known as the refractory density equation or age-structured equation:

\begin{cases} \frac{\partial}{\partial\_t} u(t,a) + \frac{\partial}{\partial\_a} u(t,a) + f(a, \xi(t)) u(t,a) = 0,\\ u(t,0) = \int\_{0}^{\infty} f(a, \xi(t)) u(t,a) da,\\ \xi(t) = \int\_{0}^{t} h(t-s) u(s,0) ds, \end{cases}

for some functions f and h. The functions I used in the script below are f(a,\xi) = 1 + \xi and h(t) = \exp(-t).

I tried `MethodOfLines.jl` and in particular this [doc](https://docs.sciml.ai/MethodOfLines/dev/tutorials/PIDE/), but did not succeed. Here is the script I tried inspired from the doc.

```julia-auto
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets

# Parameters of the PDE and initial condition
tmin = 0.0
tmax = 10.0
f(a, ξ) = 1 + ξ
h(t) = exp(t)
amin = 0.0
amax = 3.0
u_in(a) = 1.0 * (0 < a) * (a < 2)

# Parameters, variables, derivatives and integral
@parameters s t a
@variables u(..) integrand_s(..) integrand_a(..)
Dt = Differential(t)
Da = Differential(a)
Ia = Integral(a in DomainSets.ClosedInterval(amin, amax))
Is = Integral(s in DomainSets.ClosedInterval(tmin, t))

# PDE and boundary conditions
eqs = [
    Dt(u(t, a)) + Da(u(t, a)) + integrand_a(t, a) ~ 0
    integrand_s(s, t, a) ~ h(t - s) * u(t, 0)
    integrand_a(t, a) ~ f(a, Is(integrand_s(s, t, a))) * u(t, a)
]
bcs = [
    u(0, a) ~ u_in(a),
    u(t, 0) ~ Ia(integrand_a(t, a)),
    integrand_a(0, a) ~ 0,
    integrand_s(0, t, a) ~ 0
]

# Space and time domains
domains = [
    t ∈ Interval(tmin, tmax),
    s ∈ Interval(tmin, tmax),
    a ∈ Interval(amin, amax)
]

# PDE system
@named pde_system = PDESystem(eqs, bcs, domains, [t, a], [u(t, a), integrand_s(s, t, a), integrand_a(t, a)])

# Method of lines discretization
discretization = MOLFiniteDifference([a => 100, s => 100], t)
prob = discretize(pde_system, discretization)

```

At this point, the following error is raised:

```julia-auto
ERROR: ArgumentError: Integration Domain only supported for integrals from start of iterval to the variable, got 0.0 .. t in Integral(s, 0.0 .. t)(integrand_s(s, t, a))

```

If I understand correctly, it does not support the fact that the upper-bound of the integration domain is the variable t. However, there is a similar case in the doc (by the way, the doc is a bit confusing about this: in my humble opinion there are too many x's, both as integration variable and integral bound).

I recently posted this [related topic](https://discourse.julialang.org/t/solving-a-non-linear-p-i-de-neural-field-equation/132678).

---

<div class="post-metadata">

### Author: ![dawbarton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dawbarton/32/215461_2.png) [@dawbarton](https://discourse.julialang.org/u/dawbarton)
#### Post date: [September 26, 2025, 5:43pm UTC](https://discourse.julialang.org/t/solving-a-non-linear-p-i-de-age-structured-equation/132677/2 "2025-09-26T17:43:59Z")

</div>

I’m not sure that MethodOfLines.jl is entirely appropriate for this type of PDE, though it’s been a while since I’ve looked at it so I don’t know. You might want to look at the escalator boxcar train (EBT) method which solves a class of age-structured problems similar to this. IIRC it amounts to integrating along the characteristics of the PDE. I’m not aware of any Julia implementations, but it shouldn’t be too hard to do.

You have to be a little careful with some equations of this form because the hyperbolic nature of the PDE can lead to shock waves developing in the solutions. I’ve no idea about what happens in this specific case, but when I last worked on age-structured problems I got caught out by that.

---

<div class="post-metadata">

### Author: ![yolhan\_mannes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yolhan_mannes/32/220485_2.png) [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)
#### Post date: [September 26, 2025, 5:57pm UTC](https://discourse.julialang.org/t/solving-a-non-linear-p-i-de-age-structured-equation/132677/3 "2025-09-26T17:57:50Z")

</div>

Pretty sure it’s not unless you add a little bit of diffusion did you try a handmade code with maybe Integral.jl and OrdinaryDiffEq.jl around to help ? Finite differences scheme are not too hard to develop and even better a finite volume method  
I don’t have my computer but here is chatgpt take didn’t give him f I think:

```julia-auto
using OrdinaryDiffEq
using LinearAlgebra

f(a, xi) = 1 + xi
A = 10.0                
N = 400                  
Δa = A / N
a_edges = range(0, A; length = N+1)
a_centers = (a_edges[1:end-1] .+ a_edges[2:end]) ./ 2

u_in(a) = 1.0 * (0 < a) * (a < 2)
u0 = u_in.(a_centers)
ξ0 = 0.0

y0 = vcat(u0, ξ0)

function fv_rhs!(dy, y, p, t)
    u = @view y[1:N]
    xi = y[end]
    du = @view dy[1:N]
    B = sum(f(ai, xi) * ui for (ai, ui) in zip(a_centers, u)) * Δa
    uL = B
    du[1] = - (u[1] - uL) / Δa - f(a_centers[1], xi) * u[1]
    @inbounds for i in 2:N
        du[i] = - (u[i] - u[i-1]) / Δa - f(a_centers[i], xi) * u[i]
    end
    dy[end] = B - xi
    return nothing
end
tspan = (0.0, 10.0)
prob = ODEProblem(fv_rhs!, y0, tspan)
@time sol = solve(prob,BS3(); saveat=0.1,dtmin=0.00001);

```

Probably really bad 🤣 note that here it used a trick to compute the \xi in the ode

---

<div class="post-metadata">

### Author: ![yolhan\_mannes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yolhan_mannes/32/220485_2.png) [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)
#### Post date: [September 26, 2025, 8:07pm UTC](https://discourse.julialang.org/t/solving-a-non-linear-p-i-de-age-structured-equation/132677/4 "2025-09-26T20:07:23Z")

</div>

gave me  
 ![age_structured_model](https://global.discourse-cdn.com/julialang/original/3X/8/2/8220d0a324c3736616458d11ef7d5505df77e1d1.gif)

its really dependant of the number of points so I think something is hard in there

---

<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: [September 26, 2025, 9:44pm UTC](https://discourse.julialang.org/t/solving-a-non-linear-p-i-de-age-structured-equation/132677/5 "2025-09-26T21:44:38Z")

</div>

Hi Julien!

Welcome to the Julia community. In Quentin’s PhD that you know, we solved a very similar equation with a finite volume method of order 1 with upwind. If I am not mistaken, your `u` is a probability distribution and your scheme has to ensure that mass is conserved (at least). If your numerical scheme “loses” masses, it is very numerically unstable.

You can then time step using SSPRK (from [DifferentialEquations.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/DifferentialEquations)). I can try to dig the code if you need it.

---

<div class="post-metadata">

### Author: ![jucheval](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jucheval/32/213961_2.png) [@jucheval](https://discourse.julialang.org/u/jucheval)
#### Post date: [September 29, 2025, 8:51am UTC](https://discourse.julialang.org/t/solving-a-non-linear-p-i-de-age-structured-equation/132677/6 "2025-09-29T08:51:02Z")

</div>

Thank you all !  
I will have a look at your ideas.
