Solving Integro-Differential Equations with Physics-Informed Neural Networks (PINNs)

The integral of function $u(x)$,

\[\int_{0}^{t}u(x)dx\]

where $x$ is variable of integral and $t$ is variable of integro-differential equation, is defined as

using ModelingToolkit
@parameters t
@variables i(..)
Ii = Symbolics.Integral(t in DomainSets.ClosedInterval(0, t))

In multidimensional case,

Ix = Integral((x, y) in DomainSets.UnitSquare())

The UnitSquare domain ranges both x and y from 0 to 1. Similarly, a rectangular or cuboidal domain can be defined using ProductDomain of ClosedIntervals.

Ix = Integral((
    x, y) in DomainSets.ProductDomain(ClosedInterval(0, 1), ClosedInterval(0, x)))

1-dimensional example

Let's take an example of an integro-differential equation:

\[\frac{∂}{∂t} u(t) + 2u(t) + 5 \int_{0}^{t}u(x)dx = 1 \ \text{for} \ t \geq 0\]

and boundary condition

\[u(0) = 0\]

using NeuralPDE, Lux, ModelingToolkit, Optimization, OptimizationOptimJL, DomainSets
using DomainSets: Interval
using IntervalSets: leftendpoint, rightendpoint
using Plots

@parameters t
@variables i(..)
Di = Differential(t)
Ii = Integral(t in DomainSets.ClosedInterval(0, t))
eq = Di(i(t)) + 2 * i(t) + 5 * Ii(i(t)) ~ 1
bcs = [i(0.0) ~ 0.0]
domains = [t ∈ Interval(0.0, 2.0)]
chain = Lux.Chain(Lux.Dense(1, 15, Lux.σ), Lux.Dense(15, 1))

strategy_ = QuadratureTraining()
discretization = PhysicsInformedNN(chain,
    strategy_)
@named pde_system = PDESystem(eq, bcs, domains, [t], [i(t)])
prob = NeuralPDE.discretize(pde_system, discretization)
callback = function (p, l)
    println("Current loss is: $l")
    return false
end
res = Optimization.solve(prob, BFGS(); maxiters = 100)

retcode: Failure
u: ComponentVector{Float64}(layer_1 = (weight = [-0.503838251248447; 1.3612692298008706; … ; -0.1527269553167059; 1.346482665457514;;], bias = [0.10072177244335304, 0.16694398158837953, -0.32750828139645893, 1.0011958844590332, 1.3257304453655714, 0.5562265455795979, -0.25534234307580417, 0.676185362342184, 0.7601004124349487, 0.4932637225619346, -0.03774825702448208, 0.3088760674668855, -1.0004969590405883, -0.4863806633883433, -0.21279324473815744]), layer_2 = (weight = [-1.637408467282734 0.23929683599678248 … -0.07791114783256951 -2.0711283290917515], bias = [-0.9605143785696317]))

Plotting the final solution and analytical solution

ts = [leftendpoint(d.domain):0.01:rightendpoint(d.domain) for d in domains][1]
phi = discretization.phi
u_predict = [first(phi([t], res.u)) for t in ts]

analytic_sol_func(t) = 1 / 2 * (exp(-t)) * (sin(2 * t))
u_real = [analytic_sol_func(t) for t in ts]
plot(ts, u_real, label = "Analytical Solution")
plot!(ts, u_predict, label = "PINN Solution")