Non-autonomous ODE
We prove the existence of a $2\pi$-periodic solution to the non-autonomous ODE
\[\ddot{u}(t) + \beta_1 \dot{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: Formulation
using RadiiPolynomialstruct LOp{T} <: AbstractDiagonalOperator β₁ :: T β₂ :: Tendfunction 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))endF(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: Approximation (floating-point arithmetic)
The approximate zero
β_approx = [0.1, 4.0, 1.0]c = Sequence(Fourier(1, 1.0), [0.5, 0.0, 0.5])K = 10u_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 solutionNewton'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 | 5.7e+01
2 | 1.9515e-02 | 2.9687e-02 | 2.6e+01
3 | 3.2174e-03 | 3.1294e-03 | 1.6e+01
4 | 6.4577e-06 | 5.2631e-06 | 1.1e+01
5 | 2.5471e-11 | 3.6997e-11 | 8.1e+00
6 | 2.8610e-17 | 2.5831e-17 | 6.1e+00using CairoMakielines(LinRange(-π, π, 101), t -> real(u_bar(t)))
The approximate inverse
Π = interval( Projection(Fourier(K, 1.0)) )A_K = interval(inv(mid.(Π * DF(u_bar, β_approx) * Π)))A = A_K + (interval(I) - Π)Step 3: Bounds (interval arithmetic)
ν = interval(1.1) # does not have to be exactly 11/10X = Ell1(GeometricWeight(ν))β_interval = [I"0.1", I"4.0", I"1.0"]c_interval = interval(c)u_bar_interval = interval(u_bar)#- Y boundY = norm(A * F(u_bar_interval, β_interval, c_interval), X)#- Z₁ bound@assert K ≥ sup(sqrt(β_interval[2])) # |ℓ_k| is increasing beyond this pointΠ_2Kp1 = interval( Projection(Fourier(2K+1, 1.0)) )Z₁ = opnorm(Π_2Kp1 - A * DF(u_bar_interval, β_interval) * Π_2Kp1, X)#- Z₂ boundR = InfZ₂ = 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 = [8.65329e-10, 8.65333e-10]_com
│ Z₁ = [0.0115027, 0.0115028]_com
│ Z₂ = [8.94571, 8.94572]_com
│ R = Inf
│ Δ = [0.977126, 0.977127]_com
└ roots = ([8.75398e-10, 8.75402e-10]_com, [0.220999, 0.221]_com)Step 4: Conclusion
inf(ie) # smallest error8.754019290271712e-10