• Nonlinear diffusion
• • • Step 1: Problem definition
    • Step 2: Approximate zero (floating-point arithmetic)
    • Step 3: Approximate inverse (floating-point arithmetic)
    • Step 4: Estimating the bounds (interval arithmetic)

Nonlinear diffusion

In this example, we prove the existence of a steady-state of a reaction-diffusion equation with nonlinear diffusion

\[\begin{cases} \partial_t u = \Delta \Phi(u) + R(u), & x \in (0, 1), \\ \displaystyle \frac{\partial u}{\partial n} = 0, & x \in \{0, 1\}, \end{cases}\]

with

\[\begin{aligned} \Phi(u) &= u^2, \\ R(u) &= u - u^2 + g(x), \\ g(x) &= \frac{1}{2} + 3 \cos(\pi x) + 2 \cos(2 \pi x) - \cos(3 \pi x) + 6 \cos(4 \pi x). \end{aligned}\]

See the reference^[1] for more details.

Step 1: Problem definition

Φ(u) = u^2
DΦ(u) = exact(2) * u

R(u, g) = u - u^2 + g
DR(u) = exact(1) - exact(2) * u

F(u, g) = Laplacian() * Φ(u) + R(u, g)
DF(u) = Laplacian() * Multiplication(DΦ(u)) + Multiplication(DR(u))

Step 2: Approximate zero (floating-point arithmetic)

using RadiiPolynomial

g = Sequence(evensym(Fourier(4, π)), [1/2, 3/2, 1, -1/2, 3])

K = 20

u_init = Sequence(evensym(Fourier(K, π)),
    [1.362741344081890 ; 0.052107816731015 ; 0.008200891820525 ; -0.002635040629728 ; 0.007018830076427 ; zeros(K-4)])

u_bar, newton_success = newton(u -> (F(u, g), DF(u)), u_init; verbose = true)

Newton's method: Inf-norm, tol = 1.0e-12, ϵ = 2.220446049250313e-16, maxiter = 15, convergence criterion: |F(x)| ≤ tol
  Iteration       |F(x)|      |DF(x)\F(x)|      ETA (s)
-----------------------------------------------------------
      0      |  3.4102e-01  |  5.0415e-04  |    0.0e+00
      1      |  1.2425e-04  |  7.5988e-08  |    9.2e+01
      2      |  7.0232e-12  |  4.0072e-15  |    4.3e+01
      3      |  2.2204e-16  |  1.2934e-16  |    2.6e+01

using GLMakie
lines(LinRange(0, 1, 201), t -> real(u_bar(t)))

Step 3: Approximate inverse (floating-point arithmetic)

struct PseudoInverseLaplacian <: AbstractDiagonalOperator end
RadiiPolynomial.getcoefficient(::PseudoInverseLaplacian, (codom, i)::Tuple{SymmetricSpace{<:Fourier},Integer}, (dom, j)::Tuple{SymmetricSpace{<:Fourier},Integer}) =
    (i == j) & !(i == j == 0) ? inv(- (frequency(dom) * exact(i))^2) : zero(frequency(dom))

Π = Projection(evensym(Fourier(K, π)))
Π_2K = Projection(evensym(Fourier(2K, π)))
A_finite = interval(Π * inv(Π_2K * DF(u_bar) * Π_2K) * Π)

M_cos = DΦ(u_bar)
M_fou = Projection(Fourier(K, π)) * DΦ(u_bar)
M_fou_inv = inv(M_fou)
M_cos_inv_interval = interval(Π * M_fou_inv)

A_tail = (Multiplication(M_cos_inv_interval) - interval(Π) * Multiplication(M_cos_inv_interval) * interval(Π)) * PseudoInverseLaplacian()

A = A_finite + A_tail

Step 4: Estimating the bounds (interval arithmetic)

g_interval = interval(g)
u_bar_interval = interval(u_bar)

#- Y bound

Y = norm(A * F(u_bar_interval, g_interval), 1)

#- Z₁ bound

Z₁_finite = opnorm(interval(Π_2K) - A * DF(u_bar_interval) * interval(Π_2K), 1)

Z₁_tail = norm(exact(1) - M_cos_inv_interval * DΦ(u_bar_interval)) +
    norm(M_cos_inv_interval, 1) / (exact(K+1)^2 * interval(π)^2) * norm(DR(u_bar_interval), 1)

Z₁ = max(Z₁_finite, Z₁_tail)

#- Z₂ bound
opnorm_A_Delta = max(opnorm(A * Laplacian() * interval(Π), 1), norm(M_cos_inv_interval, 1))

Z₂ = exact(4) * opnorm_A_Delta

#

ie, contraction_success = interval_of_existence(Y, Z₁, Z₂, Inf; verbose = true)

┌ Info: success: interval found
│ Y = Interval{Float64}(3.5423787202270646e-12, 3.546159214468317e-12, com, true)
│ Z₁ = Interval{Float64}(0.00018710232408672816, 0.00018710232408672967, com, true)
│ Z₂ = Interval{Float64}(1.6233034498235595, 1.623303449823564, com, true)
│ R = Inf
│ Δ = Interval{Float64}(0.999625830347593, 0.9996258303476057, com, true)
└ roots = (Interval{Float64}(3.54288491950678e-12, 3.546920095952698e-12, com, true), Interval{Float64}(1.2318250143332148, 1.2318250143332226, com, true))

inf(ie) # smallest error

3.546920095952698e-12