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 ModelingToolkit, NeuralPDE, SciMLBase, Lux, Optimization, OptimizationOptimJL, DomainSets
using Optim: BFGS
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: MaxIters
u: ComponentVector{Float64}(layer_1 = (weight = [-0.4805252491153959; 1.8283248856784975; … ; -1.172235927178484; 1.2636057444234692;;], bias = [-0.8175493985542497, -0.07783687866281397, -0.5757347636670471, -0.19679231276968823, 1.0769596525897733, -0.6292596351398814, 0.5239652981427373, -0.8859280392774014, -1.481326314591075, -0.8147285208808943, -0.9748539664241878, -0.30315542390930317, -0.814798411459576, -0.7983040893417988, 0.7799307290014992]), layer_2 = (weight = [-0.5099637682557744 -0.8080411295872586 … -2.0272196813195578 1.007739792381573], bias = [0.5964119118507744]))

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")