RasmahRasmah

Physics-Informed Neural Networks

println("Julia ", VERSION)
Julia 1.12.6

Overview

A physics-informed neural network (PINN) solves a partial differential equation by training a neural network to satisfy the equation directly, without any mesh. Instead of discretizing the domain into elements and assembling a stiffness matrix, a PINN represents the solution u(x) as a small network and asks it to make the PDE residual vanish at a scattered set of collocation points.

For the Poisson problem the strong form is

\[-\\Delta u = f \\quad\\text{on }\\Omega,\\qquad u = g \\quad\\text{on }\\partial\\Omega.\]

The network is trained to minimize the interior residual Δu + f plus a soft boundary term u − g. Because the whole thing is differentiable, the training loop simply descends the residual with Rasmah’s native optimizer.

using Rasmah
using Random

A 1-D Poisson problem

The problem −u'' = π² sin(πx) on [0, 1] with u(0) = u(1) = 0 has the exact solution u(x) = sin(πx). We build a small network, then train it with pinn_solve — the source f and boundary value g are passed as plain functions.

rng = MersenneTwister(1)
net = MLP([1, 32, 32, 1])
θ0 = xavier_init(rng, net.sizes)

f(x) = π^2 * sin(π * x[1])
g(x) = 0.0

sol = pinn_solve(net, θ0, f, g, [0.0], [1.0]; n = 200, epochs = 600, lr = 1e-3, seed = 1)
PINNResult{MLP, typeof(tanh), Float64}(MLP([1, 32, 32, 1]), [-0.1412453125566969, 0.2758770822088113, -0.310677507168315, -0.4211271570989069, -0.5704008856616326, 0.36256560050245334, 0.3311452450323901, -0.25449927887966095, 0.18615583798346486, -0.2867002942124032  …  -0.1769027437980979, -0.1307153940732304, 0.4732001074890667, -0.3364298697511368, 0.24990762879638706, -0.5992500157136137, 0.12087284714860906, -0.12740930813794113, 0.1947707766885065, -0.2735110796873325], tanh, [46.51421191599097, 46.16802071494183, 45.824721785604225, 45.484119982127645, 45.145920761234336, 44.8097236260724, 44.47512557537261, 44.141756104890746, 43.80924065861861, 43.47718112766864  …  0.0005633202906389907, 0.0005576320174117083, 0.0005521334445425702, 0.0005468184516342955, 0.0005416794956048916, 0.0005367103453442835, 0.0005319057947327736, 0.0005272579734989484, 0.0005227630405631626, 0.0005184143222196258])

The result sol is callable — evaluate it anywhere. We check the maximum pointwise error against the exact solution, and print both:

err = maximum(abs(sol([x]) - sin(π * x)) for x in range(0.0, 1.0; length = 41))
println("PINN max error vs sin(πx) = ", round(err; digits = 4))
println("exact u(0.5)            = ", round(sin(π * 0.5); digits = 4))
println("PINN  u(0.5)            = ", round(sol([0.5]); digits = 4))
PINN max error vs sin(πx) = 0.0056
exact u(0.5)            = 1.0
PINN  u(0.5)            = 1.0021

Enforcing the boundary exactly

Soft boundary terms can leave a small residual at the boundary. Passing hard_bc = true instead multiplies the network by a distance factor d(x) = x(1−x), so u = d·net satisfies the Dirichlet condition exactly — the boundary loss is dropped entirely.

sol_hard = pinn_solve(net, θ0, f, g, [0.0], [1.0]; n = 200, epochs = 600, lr = 2e-3,
                      hard_bc = true, seed = 1)
println("u(0) = ", sol_hard([0.0]), "   u(1) = ", sol_hard([1.0]))
err_hard = maximum(abs(sol_hard([x]) - sin(π * x)) for x in range(0.05, 0.95; length = 19))
println("hard-BC max error = ", round(err_hard; digits = 4))
u(0) = 0.0   u(1) = 0.0
hard-BC max error = 0.0026

Helmholtz equation

The Helmholtz problem −u'' − k²u = f is solved by pinn_solve_helmholtz. With k = 1 and f(x) = (4π² − 1) sin(2πx) the exact solution is u = sin(2πx):

k = 1.0
fh(x) = (4π^2 - k^2) * sin(2π * x[1])
sol_h = pinn_solve_helmholtz(net, θ0, fh, k, g, [0.0], [1.0]; n = 200, epochs = 600, lr = 1e-3, seed = 13)
err_h = maximum(abs(sol_h([x]) - sin(2π * x)) for x in range(0.0, 1.0; length = 41))
println("Helmholtz max error vs sin(2πx) = ", round(err_h; digits = 4))
Helmholtz max error vs sin(2πx) = 0.1722

2-D elasticity with the Deep Ritz method

For vector PDEs the Deep Ritz method minimizes the variational (energy) functional rather than the strong-form residual — cheaper and better-conditioned, because it only needs first derivatives. deep_ritz_elasticity solves 2-D plane-strain elasticity with a manufactured solution u = [sin(πx)sin(πy), sin(πx)sin(πy)] and its matching body force:

λ, μ = 1.0, 1.0
ux(x) = sin(π * x[1]) * sin(π * x[2])
fx(x) = (λ + 3μ) * π^2 * sin(π * x[1]) * sin(π * x[2]) - (λ + μ) * π^2 * cos(π * x[1]) * cos(π * x[2])
f2(x) = [fx(x), fx(x)]

net2 = MLP([2, 16, 16, 2])
θ2 = xavier_init(rng, net2.sizes)
sol_e = deep_ritz_elasticity(net2, θ2, f2, λ, μ, [0.0, 0.0], [1.0, 1.0];
                             n = 200, epochs = 500, lr = 3e-3, seed = 14)

err_e = maximum(abs(sol_e([x, y])[1] - ux([x, y])) for x in range(0.1, 0.9; length = 7), y in range(0.1, 0.9; length = 7))
println("Deep Ritz elasticity max error = ", round(err_e; digits = 4))
Deep Ritz elasticity max error = 0.0824

Next steps