• Non-autonomous ODE
• • • 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)

Non-autonomous ODE

In this example, we prove the existence of a $2\pi$-periodic solution to the non-autonomous ODE

\[u''(t) + \beta_1 u'(t) + \beta_2 u(t) - u(t)^2 = \beta_3 \cos(t),\]

where $\beta_1 = 1/10$, $\beta_2 = 4$, $\beta_3 = 1$.

Step 1: Problem definition

using RadiiPolynomial

struct LOp{T} <: AbstractLinearOperator
    β₁ :: T
    β₂ :: T
end
RadiiPolynomial.domain(::LOp, codom::Fourier) = codom
RadiiPolynomial.codomain(::LOp, dom::Fourier) = dom
function RadiiPolynomial.getcoefficient(L::LOp, (codom, k)::Tuple{Fourier,Integer}, (dom, j)::Tuple{Fourier,Integer})
    x = inv(-exact(k)^2 + exact(im*k)*L.β₁ + L.β₂)
    return ifelse(k == j, x, zero(x))
end

F(u, β, c) = u + LOp(β[1], β[2]) * (-u^2 - β[3] * c)

DF(u, β) = exact(I) + LOp(β[1], β[2]) * Multiplication(-exact(2) * u)

DF (generic function with 1 method)

Step 2: Approximate zero (floating-point arithmetic)

β_approx = [0.1, 4.0, 1.0]

c = Sequence(Fourier(1, 1.0), [0.5, 0.0, 0.5])

K = 10

u_init = zeros(ComplexF64, Fourier(K, 1.0))

u_bar, success_newton = newton(u -> (F(u, β_approx, c), DF(u, β_approx)), u_init; verbose = true)
conjugacy_symmetry!(u_bar) # impose real-valued solution

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      |  1.6657e-01  |  1.6657e-01  |    0.0e+00
      1      |  1.3873e-01  |  1.3913e-01  |    4.5e+01
      2      |  1.9515e-02  |  2.9687e-02  |    2.1e+01
      3      |  3.2174e-03  |  3.1294e-03  |    1.3e+01
      4      |  6.4577e-06  |  5.2631e-06  |    8.8e+00
      5      |  2.5471e-11  |  3.6997e-11  |    6.4e+00
      6      |  2.8610e-17  |  2.5831e-17  |    4.8e+00

using GLMakie
lines(LinRange(-π, π, 101), t -> real(u_bar(t)))

Step 3: Approximate inverse (floating-point arithmetic)

Π = interval( Projection(Fourier(K, 1.0)) )

A_K = interval(inv(mid.(Π * DF(u_bar, β_approx) * Π)))

A = A_K + (interval(I) - Π)

Step 4: Estimating the bounds (interval arithmetic)

ν = interval(1.1) # does not have to be exactly 11/10
X = Ell1(GeometricWeight(ν))

β_interval = [I"0.1", I"4.0", I"1.0"]
c_interval = interval(c)
u_bar_interval = interval(u_bar)

#- Y bound

Y = norm(A * F(u_bar_interval, β_interval, c_interval), X)

#- Z₁ bound
@assert K ≥ sup(sqrt(β_interval[2]))
Π_2Kp1 = interval( Projection(Fourier(2K+1, 1.0)) )

Z₁ = opnorm(Π_2Kp1 - A * DF(u_bar_interval, β_interval) * Π_2Kp1, X)

#- Z₂ bound
R = Inf

Z₂ = exact(2) * max(opnorm(A_K, X), interval(1))

#

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

┌ Info: success: interval found
│ Y = Interval{Float64}(8.653295556856143e-10, 8.65332313699294e-10, com, true)
│ Z₁ = Interval{Float64}(0.011502747276604454, 0.011502747276604998, com, true)
│ Z₂ = Interval{Float64}(8.94571074327214, 8.945710743272164, com, true)
│ R = Inf
│ Δ = Interval{Float64}(0.977126803159674, 0.977126803159725, com, true)
└ roots = (Interval{Float64}(8.753989628760099e-10, 8.754019290271712e-10, com, true), Interval{Float64}(0.22099915304131024, 0.2209991530413138, com, true))

inf(ie) # smallest error

8.754019290271712e-10