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

Specific Domains Numerics
,
1

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, but did not succeed. Here is the script I tried inspired from the doc.

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:

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.

2

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.

3

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:

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 :rofl: note that here it used a trick to compute the \xi in the ode

1 Like
4

gave me
age_structured_model

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

5

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). I can try to dig the code if you need it.

6

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