How to use a convolution integral with NeuralPDE.discretize

Hello everyone!

I am trying to include the convolution integral

\int_0^tKr(t -\tau)\dot x(\tau)d\tau

in my equation and I keep coming up with this error when trying to do NeuralPDE.discretize

ERROR: KeyError: key :tau not found

Below is a reduced version of my code.

using NeuralPDE, ModelingToolkit, Optimization, OptimizationOptimJL
using DomainSets
using Lux: Chain, Dense

K_func(t) = exp(-t)*cos(t)

@parameters t 
@parameters tau
@variables x(..)

Dt = Differential(t)
DDt = Differential(t)^2

integral_tau = Integral(tau in DomainSets.ClosedInterval(0.0, t))

eq = DDt(x(t)) + integral_tau(K_func(t - tau) * Dt(x(tau))) + x(t) ~ 0

domains = [
    t ∈ Interval(0.0, 60.0),
    tau ∈ Interval(0.0, 60.0)
]

bcs = [
    x(0.0) ~ 0.0,
    Dt(x(0.0)) ~ 0.0,
]

chain = Chain([Dense(1, 4), Dense(4, 1)])
strategy = QuasiRandomTraining(10)
discretization = PhysicsInformedNN(chain, strategy)
@named pde_system = PDESystem(eq, bcs, domains, [t], [x])

prob = NeuralPDE.discretize(pde_system, discretization) # KeyError: key :tau not found

opt = Optim.LBFGS()
maxiters=1
res = Optimization.solve(prob, opt; maxiters)

This happens with julia v1.10.11 and the latest version right now, v1.12.6.
This is on NeuralPDE v6.0.0

I’ve tried both @parameters and @variables; I understand that @parameters for independent variables and @variables are for functions of variables. I’m not sure where tau would fall under the current implementation.

I got a version running when I include tau as one of the independent variables, but that required also rewriting every instance of x(t) to also include tau “x(t, tau)”.

There was a somewhat similar issue a few years ago (I couldn’t find how to do the linked topic sorry, change the underscore to a dot, https://discourse.julialang_org/t/way-to-write-the-integral-with-modelingtoolkit/82347/5), so I am hoping that this is an error on my part.

Any help would be greatly appreciated, thanks in advance!

Similar issue from a few years ago, but doesn’t seem to have a resolution either.

How WOULD your typical ODE solver be able to handle this type of integral, if the integral couldn’t some how be re-written as a system of ODEs? Because the step the integrator will take is not just dependent on the current state, but also the history of the state — the x(tau) term where tau ranges from time 0 to current time t.

But any ODE interface, you write a function that takes the current state x and returns the rates of change x_dot. You don’t really have access to the state history, and although I suppose you could hack it in via some global variables
your ODE integrator expects, AFAIK, for those rates to be dependent ONLY on the current state.

@mark.garnett I understand the issue you’re pointing out. Looking more into this landed me onto a MATLAB thread that describes a similar equation.

In the end they solved it in the ways you described. One of the latest comment chains describes how to rewrite as a coupled ODE under certain conditions. Another comment chain includes a few more solutions; namely a Laplace transform or storing the history.

In cases where the ODE can’t be rewritten; is it up to the user to store the history (using global variables or otherwise) like the MATLAB problem or is there an implementation I am missing? This also brings a different issue to mind, how does/could NeuralPDE handle this? To my understanding it doesn’t take steps in the same way a typical ODE solver does while training the Lux model.

I don’t know what to do if you can’t re-write the integral-differential equation as an ODE in the context of Modeling Toolkit. You don’t have access to the history of the variables AFAICT, maybe some of MTK devs have an idea of how.

I ran into this situation where the problem was heat diffusion into a semi-infinite space, where the surface temperature was a function of time, and because there was a 1/sqrt(t-tau) term in the integral that was beyond my mathematical powers to handle turn into a system of ODEs.

It seems like you do. Here is an example from NeuralPDE of solving an integro differential equation.

Great! Thank you. I’ll see if I can make it work for the problem I was trying to solve.

Briefly reading through NeuralPDE’s paper, the section describing Integro-Differential Equations gives an answer to how NeuralPDE.jl handles the history problem.

“Because the neural network defines the value of 𝑢(𝑥) over all space, one can simply use a quadrature technique, such as those in Quadrature.jl , to evaluate the current value of the integral and thus utilize these terms within the PDE equations.”

While the paper is outdated, the example that @cmichelenstrofer gave is from the up-to-date docs. In that case this should be a syntax issue then, @mark.garnett let us know if you have better luck than I did when writing the equation for your problem!