Nonlinear elliptic system of PDEs

We can also solve nonlinear systems, such as the system of nonlinear elliptic PDEs

\[\begin{aligned} \frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2} &= uf(\frac{u}{w}) + \frac{u}{w}h(\frac{u}{w}) \\ \frac{\partial^2w}{\partial x^2} + \frac{\partial^2w}{\partial y^2} &= wg(\frac{u}{w}) + h(\frac{u}{w}) \\ \end{aligned}\]

where f, g, h are arbitrary functions. With initial and boundary conditions:

\[\begin{aligned} u(0,y) &= y + 1 \\ w(1,y) &= [\cosh(\sqrt[]{f(k)}) + \sinh(\sqrt[]{f(k)})]\cdot(y + 1) \\ w(x,0) &= \cosh(\sqrt[]{f(k)}) + \sinh(\sqrt[]{f(k)}) \\ w(0,y) &= k(y + 1) \\ u(1,y) &= k[\cosh(\sqrt[]{f(k)}) + \sinh(\sqrt[]{f(k)})]\cdot(y + 1) \\ u(x,0) &= k[\cosh(\sqrt[]{f(k)}) + \sinh(\sqrt[]{f(k)})] \\ \end{aligned}\]

where k is a root of the algebraic (transcendental) equation f(k) = g(k).

This is done using a derivative neural network approximation.

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

@parameters x, y
Dx = Differential(x)
Dy = Differential(y)
@variables Dxu(..), Dyu(..), Dxw(..), Dyw(..)
@variables u(..), w(..)

# Arbitrary functions
f(x) = sin(x)
g(x) = cos(x)
h(x) = x
root(x) = f(x) - g(x)

# Analytic solution
k = find_zero(root, (0, 1), Bisection())                            # k is a root of the algebraic (transcendental) equation f(x) = g(x)
θ(x, y) = (cosh(sqrt(f(k)) * x) + sinh(sqrt(f(k)) * x)) * (y + 1)   # Analytical solution to Helmholtz equation
w_analytic(x, y) = θ(x, y) - h(k) / f(k)
u_analytic(x, y) = k * w_analytic(x, y)

# Nonlinear Steady-State Systems of Two Reaction-Diffusion Equations with 3 arbitrary function f, g, h
eqs_ = [
    Dx(Dxu(x, y)) + Dy(Dyu(x, y)) ~
    u(x, y) * f(u(x, y) / w(x, y)) +
    u(x, y) / w(x, y) * h(u(x, y) / w(x, y)),
    Dx(Dxw(x, y)) + Dy(Dyw(x, y)) ~ w(x, y) * g(u(x, y) / w(x, y)) + h(u(x, y) / w(x, y))]

# Boundary conditions
bcs_ = [u(0, y) ~ u_analytic(0, y),
    u(1, y) ~ u_analytic(1, y),
    u(x, 0) ~ u_analytic(x, 0),
    w(0, y) ~ w_analytic(0, y),
    w(1, y) ~ w_analytic(1, y),
    w(x, 0) ~ w_analytic(x, 0)]

der_ = [Dy(u(x, y)) ~ Dyu(x, y),
    Dy(w(x, y)) ~ Dyw(x, y),
    Dx(u(x, y)) ~ Dxu(x, y),
    Dx(w(x, y)) ~ Dxw(x, y)]

bcs__ = [bcs_; der_]

# Space and time domains
domains = [x ∈ Interval(0.0, 1.0), y ∈ Interval(0.0, 1.0)]

# Neural network
input_ = length(domains)
n = 15
chain = [Chain(Dense(input_, n, σ), Dense(n, n, σ), Dense(n, 1)) for _ in 1:6] # 1:number of @variables

strategy = GridTraining(0.01)
discretization = PhysicsInformedNN(chain, strategy)

vars = [u(x, y), w(x, y), Dxu(x, y), Dyu(x, y), Dxw(x, y), Dyw(x, y)]
@named pdesystem = PDESystem(eqs_, bcs__, domains, [x, y], vars)
prob = NeuralPDE.discretize(pdesystem, discretization)
sym_prob = NeuralPDE.symbolic_discretize(pdesystem, discretization)

pde_inner_loss_functions = sym_prob.loss_functions.pde_loss_functions
bcs_inner_loss_functions = sym_prob.loss_functions.bc_loss_functions[1:6]
approx_derivative_loss_functions = sym_prob.loss_functions.bc_loss_functions[7:end]

callback = function (p, l)
    if p.iter % 10 == 0
        println("loss: ", l)
        println("pde_losses: ", map(l_ -> l_(p.u), pde_inner_loss_functions))
        println("bcs_losses: ", map(l_ -> l_(p.u), bcs_inner_loss_functions))
        println("der_losses: ", map(l_ -> l_(p.u), approx_derivative_loss_functions))
    end
    return false
end

res = solve(prob, BFGS(); maxiters = 100, callback)

phi = discretization.phi

# Analysis
xs, ys = [leftendpoint(d.domain):0.01:rightendpoint(d.domain) for d in domains]
depvars = [:u, :w]
minimizers_ = [res.u.depvar[depvars[i]] for i in 1:2]

analytic_sol_func(x, y) = [u_analytic(x, y), w_analytic(x, y)]
u_real = [[analytic_sol_func(x, y)[i] for x in xs for y in ys] for i in 1:2]
u_predict = [[phi[i]([x, y], minimizers_[i])[1] for x in xs for y in ys] for i in 1:2]
diff_u = [abs.(u_real[i] .- u_predict[i]) for i in 1:2]
ps = []
for i in 1:2
    p1 = plot(xs, ys, u_real[i], linetype = :contourf, title = "u$i, analytic")
    p2 = plot(xs, ys, u_predict[i], linetype = :contourf, title = "predict")
    p3 = plot(xs, ys, diff_u[i], linetype = :contourf, title = "error")
    push!(ps, plot(p1, p2, p3))
end

loss: 3.7112024918111697
pde_losses: [0.14104410042162774, 0.1436451181930165]
bcs_losses: [0.1730106776003898, 1.2749202218627107, 0.20700224928726008, 0.2882179279144713, 0.6505121291511782, 0.35774112986298084]
der_losses: [0.17932522841211704, 0.10462443866035225, 0.04771238813678825, 0.14344688230827693]
loss: 1.501054746739348
pde_losses: [0.29474264440741754, 0.28823326309068886]
bcs_losses: [0.0020344341816237767, 0.2002509373592204, 0.1219649322360673, 0.03312861414656543, 0.1502052354258786, 0.08988116239996251]
der_losses: [0.046912140520698364, 0.06913399517907486, 0.1113752741686876, 0.0931921136234629]
loss: 0.8751942106904196
pde_losses: [0.10351056481199512, 0.17832664420519673]
bcs_losses: [0.027980967003521116, 0.12489988662825258, 0.026838567691742633, 0.0075550919892798515, 0.16147384045483312, 0.019108895350721827]
der_losses: [0.06080530487153342, 0.04934244350756015, 0.054735409213667124, 0.06061659496211588]
loss: 0.32598825112112023
pde_losses: [0.03895951765513956, 0.027001369537130792]
bcs_losses: [0.0036550198472056314, 0.05214396922778901, 0.026122016568153904, 0.001820621421719065, 0.03129470929558137, 0.018475946611776345]
der_losses: [0.027642411038616727, 0.06364605149166935, 0.018998168446753907, 0.01622844997958457]
loss: 0.16665600286903187
pde_losses: [0.018882071279506817, 0.012687181402563178]
bcs_losses: [0.018937947918689843, 0.0344126172931111, 0.011743649648265655, 0.02882977488977046, 0.0024294767480382267, 0.00506360692270232]
der_losses: [0.008742290208631882, 0.009861598183362895, 0.0032902709917421217, 0.011775517382647361]
loss: 0.0862686470682999
pde_losses: [0.007478818138988957, 0.0032474133585276553]
bcs_losses: [0.005192537188339587, 0.009190504497947649, 0.009357899865724184, 0.0026150587012216805, 0.008721299261311545, 0.003250035667621551]
der_losses: [0.0037747173496150994, 0.020899048682114698, 0.0031584516760891092, 0.009382862680798193]
loss: 0.057413752133431396
pde_losses: [0.0035462629333717654, 0.003240186025465884]
bcs_losses: [0.0035667696859627337, 0.004368211817213553, 0.005999044995080804, 0.005292112415774677, 0.005433525002444647, 0.0031756986303185853]
der_losses: [0.0015734287997882723, 0.0056873286234851664, 0.002780193646183897, 0.012750989558341414]
loss: 0.03948062655339962
pde_losses: [0.00148630151982839, 0.0021893631518194352]
bcs_losses: [0.003216796316669018, 0.001048796338910346, 0.0026773802216137656, 0.005330266017564984, 0.0041601574130230735, 0.004781815139738422]
der_losses: [0.0008470627080889564, 0.0038913880641088575, 0.0036050477606650263, 0.006246251901369351]
loss: 0.02983748012063122
pde_losses: [0.0010024226808791033, 0.004599836042835068]
bcs_losses: [0.002074164655446723, 0.0008742854969735121, 0.0018939211737148223, 0.003733326701235049, 0.001968129915882045, 0.004404953623974432]
der_losses: [0.0010409472498695124, 0.0020772272664989107, 0.0032292466149251705, 0.002939018698396874]
loss: 0.016350532901906717
pde_losses: [0.00033692086501973053, 0.0021615583117246453]
bcs_losses: [0.001118081686297608, 0.00044069368385752793, 0.0019580231208759314, 0.0015410896212295717, 0.0004109533755493097, 0.0024748956743001044]
der_losses: [0.0008814149104771575, 0.0022548012344491804, 0.0017151131049401116, 0.001056987313185838]

ps[1]

ps[2]