Solving coupled system of linear Volterra integral equations of the second kind

[timmyfaraday]

Hi,

I'm trying to solve a coupled system of linear Volterra integral equations of the second kind h1(t) and h2(t) as part of solving a semi-Markov problem, to determine state probabilities p1(t) and p2(t). In this light, I've tried to use Integrals.jl, but I'm running into a StackOverflowError.

Is there any fix to avoid this? Or is the package Integrals.jl not suited for solving this type of problems? Below you may find a 'working' example.

Thanks in advance for your time and effort!

using Distributions, Integrals

# state 1
D₁₂     = Exponential(2.0)
f₁₂(φ)  = pdf(D₁₂, φ)
F₁(φ)   = cdf(D₁₂, φ)

# state 2
D₂₁     = Exponential(1.0)
f₂₁(φ)  = pdf(D₂₁, φ)
F₂(φ)   = cdf(D₂₁, φ)

# transition frequency densities
h₁(t) = solve(IntegralProblem((τ, p) -> f₂₁(t-τ) .* h₂.(τ), 0.0, t), QuadGKJL()).u
h₂(t) = solve(IntegralProblem((τ, p) -> f₁₂(t-τ) .* h₁.(τ), 0.0, t), QuadGKJL()).u .+ f₁₂.(t)

# state probabilities
p₁(t) = solve(IntegralProblem((τ, p) -> h₁.(τ) * F₁(t-τ), 0.0, t), QuadGKJL()).u .+ F₁(t) 
p₂(t) = solve(IntegralProblem((τ, p) -> h₂.(τ) * F₂(t-τ), 0.0, t), QuadGKJL()).u

p₁(1.0)

[lxvm]

Hi,

The coupled equations you wrote for h1 and h2 are defined recursively, which is why you have the stack overflow. Integrals.jl currently works only for integrands you can evaluate directly. Maybe some of the other SciML libraries can be useful if the integral equation can be rewritten as a differential equation. For example, if the coupled equations for h1 and h2 could be solved with an ODE solver, then p1 could integrate that solution with Integrals.jl

[ChrisRackauckas]

This kind of equation needs to use a Delay Differential Equation solver. A purely explicit integral solver cannot solve such equations. Check the DDE solver documentation for details on how to use it. In that solver, h(p,t) is a fully continuous object is continuous object with all of the past history, and you can use Integrals.jl inside of the DDE.